SR_SB_Veg_Analysis
About
Seed rain, seed bank, and vegetational dynamics of north-central Missouri remnant and reconstructed tallgrass prairies
Background:
Author: Katherine Wynne (Wynnekat@msu.edu)
Co-authors: E. Eyerly, L. Sullivan
Created: 14 June 2023
Files:
Cleaned_Spring_Cover_Data_July_2019.xlsx - contains spring p/a cover data
Cleaned_Fall_Cover_Data_September_2019.xlsx - contains fall % cover data
Seed_Rain_Data_2019_FINAL.xlsx - contains seed rain data
Seed_Bank_2020_Cleaned.xlsx - contains seed bank data
SR_Traits_Dataset.xlsx - contains trait data (work in progress)
R Version: R 4.2.2
RStudio Version: 2023.03.0+386
Package Version:
Last updated: 16 June 2023
—–
Setup
Libraries
### --- Dataset importing and exporting -----
# Import in excel files
library(readxl)
# Export dataframes to excel files
library(writexl)
### ---- Data manipulation and visualization ----
# General data cleaning, manipulation, and visualization
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.3 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors#
library(vegan)## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4library(goeveg)## This is GoeVeg 0.6.5 - build: 2023-06-06#
library(maditr)##
## To select columns from data: columns(mtcars, mpg, vs:carb)
##
##
## Attaching package: 'maditr'
##
## The following objects are masked from 'package:dplyr':
##
## between, coalesce, first, last
##
## The following object is masked from 'package:purrr':
##
## transpose
##
## The following object is masked from 'package:readr':
##
## cols#
library(cowplot)##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:lubridate':
##
## stamp#
library(BiodiversityR)## Loading required package: tcltk
## BiodiversityR 2.15-3: Use command BiodiversityRGUI() to launch the Graphical User Interface;
## to see changes use BiodiversityRGUI(changeLog=TRUE, backward.compatibility.messages=TRUE)#
library(ggrepel)
#
library(vegan3d)
### ---- Model fitting and analysis ----
# Mixed-models
## Fits mixed models
library(lme4)## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack ## Functions to analyze significance and R^2 components of mixed-models
library(lmerTest)##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step ##
library(car)## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some# Summary stats
##
library(mosaic)## Registered S3 method overwritten by 'mosaic':
## method from
## fortify.SpatialPolygonsDataFrame ggplot2
##
## The 'mosaic' package masks several functions from core packages in order to add
## additional features. The original behavior of these functions should not be affected by this.
##
## Attaching package: 'mosaic'
##
## The following objects are masked from 'package:car':
##
## deltaMethod, logit
##
## The following object is masked from 'package:lmerTest':
##
## rand
##
## The following object is masked from 'package:lme4':
##
## factorize
##
## The following object is masked from 'package:Matrix':
##
## mean
##
## The following object is masked from 'package:cowplot':
##
## theme_map
##
## The following objects are masked from 'package:goeveg':
##
## deg2rad, rad2deg
##
## The following object is masked from 'package:permute':
##
## shuffle
##
## The following objects are masked from 'package:dplyr':
##
## count, do, tally
##
## The following object is masked from 'package:purrr':
##
## cross
##
## The following object is masked from 'package:ggplot2':
##
## stat
##
## The following objects are masked from 'package:stats':
##
## binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
## quantile, sd, t.test, var
##
## The following objects are masked from 'package:base':
##
## max, mean, min, prod, range, sample, sum##
library(MASS)##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## selectlibrary(MuMIn)
library(performance)
library(see)
library(patchwork)##
## Attaching package: 'patchwork'
##
## The following object is masked from 'package:MASS':
##
## area
##
## The following object is masked from 'package:cowplot':
##
## align_plotsImport Data
spring_cover <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Vegetative Cover Data 2019/Cleaned_Spring_Cover_Data_July_2019.xlsx")
fall_cover <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Vegetative Cover Data 2019/Cleaned_Fall_Cover_Data_September_2019.xlsx")
seed_rain <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Seed Rain 2019/Seed_Rain_Data_2019_FINAL.xlsx")
seed_bank <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Seed Bank 2020/Seed_Bank_2020_Cleaned.xlsx")
traits <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Traits and Species Info/SR_SB_VEG_Traits_Dataset.xlsx")—–
Methods
Data Collection & Experimental Design
Study Sites
Our study sites were at Tucker Prairie Research Area (38○56’53.6” N, 91○59’40.0” W, Callaway County, MO) and at Prairie Fork Conservation Area (38○58’29.7” N, 91○44’03.3” W, Callaway County, MO). Tucker Prairie is a 59-hectare tract of unplowed North American tallgrass claypan prairie. Less than 0.5 % of intact tallgrass prairie ecosystem (i.e., never-been-plowed) remains in Missouri, and Tucker Prairie represents the last sizable remnant prairie in north central Missouri (Samson & Knopf, 1994). More than 250 species of plants inhabit Tucker Prairie, with representatives from 57 families and over 150 genera (Tropicos, n.d.). Native C4 grasses dominate the landscape including Sorghastrum nutans, Andropogon gerardii, Schizachyrium scoparium, and Sporobolus heterolepis (Drew, 1947; Rabinowitz & Rapp, 1980). Before 1957, Tucker Prairie was used for cattle grazing and occasional haying (Drew, 1947). From 1958 to 2002, Tucker Prairie was burned once every four years in the late winter or early spring (Rabinowitz & Rapp, 1980). Since 2002, Tucker Prairie has been managed on a 5-year burn rotation, where units are burned once in the late winter to early spring (Jan. – Mar.) and again 2-3 years later in the late summer to early fall (Aug. – Oct.).
Prairie Fork Conservation Area (PFCA) is over 450 hectares of former agricultural land being restored to tallgrass prairie and savanna ecosystems (Navarrete-Tindall et al., 2007; Newbold et al., 2019). From 2004 onwards, 16-25 hectares are newly seeded each year with native prairie species collected from Tucker Prairie and other nearby remnant prairies (Newbold et al., 2019). As a result, the plots within PFCA represent a chronosequence of reconstructed tallgrass prairies that are mainly comprised of Tucker Prairie descendants. Prior to seeding, glyphosate-resistant crops (i.e., corn and soybeans) are planted and harvested for at least three years to reduce the existing weed seed bank (Newbold et al., 2019). A seed mix containing an average of 179 species of native forbs, legumes, grasses, rushes, and sedges are broadcast seeded at a rate of 13.4 to 18.2 kg/ha once during the dormant season (Newbold et al., 2019). All mixes are comprised of at least 75% of the same species each year that are representative of nearby remnant prairies such as Tucker Prairie (Newbold et al., 2019). Reconstructions are managed using a 2-4-year burn schedule and annual spraying of invasive species (Newbold et al., 2019). For additional details on PFCA management, see Newbold et al. 2019. To capture changes in seed rain dynamics during the restoration process, we grouped our reconstructed sites into three categories, as defined by Newbold et al. (2019) as being compositionally different. We measured the seed rain in an old reconstruction (seeded in 2004), middle aged reconstruction (seeded in 2012 and 2013), and a young reconstruction (seeded in 2017).
Seed rain sampling
In May 2019, we deployed artificial turf grass seed traps (0.1 x 0.1 m) in Tucker Prairie and at each of the focal PFCA restorations. We used artificial turf grass traps instead of sticky traps due to their durability and resistance to freezing (Molau & Per Mølgaard, 1996; Bass, 2015). Turfgrass traps also do not lose effectiveness over time, like sticky traps, which are prone to contamination by soil, non-seed plant debris, and dead insects, making seed recovery difficult (Arruda et al., 2020). Furthermore, unlike funnel traps, they are easy to install and collect with minimal disturbance to surrounding vegetation and soil and do not fill with water (Schott, 1995). At each site, we randomly established ten, 5 m long transects, placed a minimum of 50 m apart from each other. We then placed traps at 1 m intervals along each transect and affixed them to the ground using ground staples (n = 5 traps per transect, 50 traps per site). We collected and replaced seed traps every 2-weeks during the growing season from May 31st to December 12th, 2019. Hereafter we refer to all captured diaspores as “seeds” despite catching a variety of fruit types (e.g., achene). After collection, all traps from a transect and sampling period were grouped together and sieved through a series of soil sieves (1 mm, 500 mm, 250 mm mesh). From these layers, we then picked out and identified seeds to the lowest taxonomic level possible using a reference collection of species inhabiting our study areas and identification guides (Coons et al., n.d.; A. C. Martin & Barkley, 1961).
Because of their extremely small size (0.3 – 0.5 mm; Yatskievych, 1999), we estimated the number of rush seeds (Juncus sp.) when there were over 200 seeds present in a sample. We first sieved the samples through 180 mm and 150 mm mesh soil sieves. Then we subsampled each sieved layer by calculating the average number of rush seeds per 1 cm2 area (n = 3) and multiplying that average by the total area covered by the sample layer. Lastly, we summed the number of estimated rush seeds per layer to calculate the total number of rush seeds in a sample.
Vegetation sampling
In 2019, at the same transects, we sampled species cover twice, once in the early summer (June - July) and again at peak biomass (September). During the first sampling period, we only recorded species presence. At peak biomass, we also sampled percent aerial species cover in a 1m2 area centered around each seed rain trap. We identified species according to Yatskievych (1999).
Seed bank sampling
In March 2020, before the growing season, we collected soil samples to assess the soil seed bank at each focal prairie. To compare the soil seed bank and the previous year’s seed rain, we collected 5 (10 x 10 x 10 cm3; 1000 cm3) soil cores approximately 0.5 m away from where we captured seed rain at each transect in 2019. We allowed the soil samples to dry before removing all non-seed plant material (e.g., roots and rhizomes) to ensure seedlings only germinated from seeds. We then homogenized all samples collected from a particular transect and subsampled ~1500 cm3 of soil to spread ~1 cm deep over sterile potting soil in plastic germination trays. Starting July 2020, we placed the trays (n = 40) in a greenhouse and watered them when dry. Amongst the 40 trays containing soil samples, we also placed an additional 10 trays containing only sterile potting soil to assess whether any contamination occurred from external sources. Once identifiable, we removed seedlings from trays. After germination ceased in July 2021, we placed all trays into a cold room for ~4 months to replicate conditions necessary to break dormancy for any still dormant seeds. Post-vernalization, we returned the trays to the greenhouse, stirred the soil samples, and monitored for any additional germination. We ended the study one year post-vernalization.
—–
General Dataset Cleaning
# Make the six letter species codes for the spring and fall cover uppercase
spring_cover[[3]] <- toupper(spring_cover[[3]])
fall_cover[[4]] <- toupper(fall_cover[[4]])
# Remove unnecessary columns from datasets
### Spring cover - Only keep Site, Transect, SPP6, Scientific Name
spring_cover <- spring_cover[,-c(5,6)]
### Fall cover - Keep everything
### Seed rain - Only keep Plot, Transect, Sampling Date, Week, SPP6, Number, Unknown
seed_rain <- seed_rain[, -c(7,8)]
### Seed bank - Keep everything
# Change names of columns to better reflect datasets (e.g. seed rain: number -> Number_Seeds)
colnames(seed_bank) <- c("Site", "Transect", "SPP6", "Number_Seedlings")
colnames(seed_rain) <- c("Site", "Transect", "Sampling_Date", "Week", "SPP6", "Number_Seeds", "Unknown")
# Change transect PFCA 2 - 10A to just 10 for the cover datasets
spring_cover <- spring_cover %>%
mutate_if(is.character,
str_replace_all, pattern = "10A", replacement = "10")
fall_cover <- fall_cover %>%
mutate_if(is.character,
str_replace_all, pattern = "10a", replacement = "10")
# Make Transect and Plot variables factors for all datasets
## Spring cover
spring_cover$Transect <- as.factor(spring_cover$Transect)
## Fall cover
fall_cover$Transect <- as.factor(fall_cover$Transect)
fall_cover$Plot <- as.factor(fall_cover$Plot)
## Seed Bank
seed_bank$Transect <- as.factor(seed_bank$Transect)
## Seed Rain
seed_rain$Transect <- as.factor(seed_rain$Transect)Lump everything to transect level data (total number of seeds per transect, total number of seedlings per transect, total amount of fall cover per transect)
####### Clean Spring
# 0.) Make a frequency column
spring_cover$Cover <- rep(1, nrow(spring_cover)) ####### Clean Fall
# 0.) Make anything less than 1 = 1
fall_cover <- fall_cover %>%
mutate(Cover = ifelse(Cover < 1, 1, Cover))
# 1.) get a sum per species per transect of cover.
fall_sum <- fall_cover %>%
group_by(Site, Transect, SPP6, Scientific_Name) %>%
summarize(Cover = sum(Cover))## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.fall_sum## # A tibble: 1,215 × 5
## # Groups: Site, Transect, SPP6 [1,214]
## Site Transect SPP6 Scientific_Name Cover
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 ACAVIR Acalypha virginica 11.5
## 2 PFCA 1 1 AMBART Ambrosia artemisiifolia 6
## 3 PFCA 1 1 BARVUL Barbarea vulgaris 3.5
## 4 PFCA 1 1 CHAFAS Chamaecrista fasciculata 388
## 5 PFCA 1 1 CONCAN Conyza canadensis 1
## 6 PFCA 1 1 ERISPP Erigeron spp. 2
## 7 PFCA 1 1 ERIVIL Eriochloa villosa 2
## 8 PFCA 1 1 EUPSER Eupatorium serotinum 3
## 9 PFCA 1 1 KUMSTI Kummerowia stipulacea 6
## 10 PFCA 1 1 LESCAP Lespedeza capitata 1
## # ℹ 1,205 more rows####### Seed Rain
# 1.) get a sum per species per transect of cover.
seed_rain_sum <- seed_rain %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seeds = sum(Number_Seeds))## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.# Clean in the same way you did the seed rain analysis
for(i in 1:nrow(seed_rain_sum)){
if(seed_rain_sum[i,3] == "AGAFAS"){seed_rain_sum[i,3] <- "AGASPP"}
if(seed_rain_sum[i,3] == "EUPSER"){seed_rain_sum[i,3] <- "EUPSPP"}
if(seed_rain_sum[i,3] == "GENPUB"){seed_rain_sum[i,3] <- "GENSPP"}
}####### Seed Bank
# 1.) get a sum per species per transect of cover.
seed_bank_sum <- seed_bank %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seedlings = sum(Number_Seedlings))## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.Remove Unknowns from Data sets
####### Remove Spring Unknowns
spring_cover <- filter(spring_cover, !grepl('UNK', SPP6))
# Looks good
list(unique(spring_cover$SPP6))## [[1]]
## [1] "ACAVIR" "ACERUB" "ACHMIL" "AGRGIG" "AGRHYE" "ALOCAR" "AMBART" "AMBPSI"
## [9] "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB" "BAPAUS" "BAPBRA"
## [17] "BARVUL" "BIDARI" "BLECIL" "BROJAP" "CAMRAD" "CAPBUR" "CARSPP" "CARBIC"
## [25] "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN" "CORLAN" "CORPAL" "CORSPP"
## [33] "CORTRI" "CROSAG" "CYPACU" "CYPECH" "DALPUR" "DESILL" "DESSES" "DESSPP"
## [41] "DICLAN" "DIGISC" "ECHPAL" "ELETEN" "ELYVIR" "ERISPP" "ERYYUC" "EUPCOR"
## [49] "EUPSER" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA" "GALOBT" "GENPUB"
## [57] "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP" "HYPPUN" "JUNBRA" "JUNMAR"
## [65] "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER" "LEPVIR" "LESCAP" "LESCUN"
## [73] "LESVIR" "LIAPYC" "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP"
## [81] "MELSPP" "MONFIS" "MUHSPP" "MYOVER" "OENBIE" "OXADIL" "PARINT" "PENDIG"
## [89] "PLAVIR" "POAPRA" "POLSAN" "POLVER" "POTSIM" "PYCPIL" "PYCTEN" "RANABO"
## [97] "RATPIN" "RHUCOP" "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI"
## [105] "SALAZU" "SCHSCO" "SENMAR" "SETFAB" "SETPAR" "SILANT" "SILINT" "SILLAC"
## [113] "SISSPP" "SOLALT" "SOLCAR" "SOLNEM" "SOLRIG" "SORNUT" "SPHOBT" "SPOHET"
## [121] "STRLEI" "SYMERI" "SYMLAN" "SYMLAT" "SYMNOV" "SYMOBL" "SYMOOL" "SYMORB"
## [129] "SYMPRA" "SYMSPP" "THLARV" "OENSPP" "TRAOHI" "TRIPER" "TRIREP" "TRISPP"
## [137] "ULMSPP" "SCISPP" "VERARV" "VERHEL" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT"
## [145] "ZIZAUR" "SCLTRI"####### Remove Fall Unknowns
fall_sum <- filter(fall_sum, !grepl('UNK', SPP6))
# Looks good
list(unique(fall_sum$SPP6))## [[1]]
## [1] "ACAVIR" "AMBART" "BARVUL" "CHAFAS" "CONCAN" "ERISPP" "ERIVIL" "EUPSER"
## [9] "KUMSTI" "LESCAP" "PENDIG" "PLAVIR" "RUDHIR" "SETFAB" "SOLALT" "STRLEI"
## [17] "SYMPIL" "ACHMIL" "CORTRI" "DALPUR" "DIGISC" "EUTGYM" "JUNSPP" "PYCTEN"
## [25] "RATPIN" "RUDSUB" "SCISPP" "SOLSPE" "SORNUT" "SYMERI" "SYMNOV" "SYMORB"
## [33] "TRICAM" "BROJAP" "HELMOL" "HYPPUN" "OENBIE" "SOLRIG" "SYMLAN" "AGRHYE"
## [41] "KUMSTR" "LESCUN" "LIAPYC" "PARINT" "SOLCAR" "TAROFF" "BAPALB" "ERASPE"
## [49] "JUNVIR" "LESVIR" "SYMPRA" "VERARV" "CERSPP" "ECHCRU" "OXADIL" "PLAMAJ"
## [57] "SILINT" "SILLAC" "VERSPP" "CORLAN" "CYPECH" "CYPSTR" "EUPCOR" "MEDLUP"
## [65] "SCHSCO" "BIDARI" "FESPAR" "MONFIS" "RANABO" "SOLNEM" "OENFIL" "ULMPUM"
## [73] "DESSPP" "ELYCAN" "AGATEN" "BAPBRA" "CARSPP" "CIRALT" "ERYYUC" "GENAND"
## [81] "KOEMAC" "LACSER" "PYCPIL" "RUMCRI" "SPOHET" "SYMOBL" "VERHEL" "BLECIL"
## [89] "BOLAST" "CYPACU" "SYMLAT" "SYMOOL" "ECHPAL" "MELSPP" "SYMANO" "SYMLAE"
## [97] "CORPAL" "ANDGER" "EUPALT" "POAPRA" "ALOCAR" "LIASPP" "LYTALA" "SPOCOM"
## [105] "SPILAC" "PANVIR" "SENMAR" "TRIFLA" "HELSPP" "ELYVIR" "HELHEL" "MYOVER"
## [113] "SALAZU" "TRAOHI" "BAPAUS" "LESPRO" "DESSES" "ROSCAR" "AMOCAN" "LEPDEN"
## [121] "CAMRAD" "RUBSPP" "ACERUB" "ZIZAUR" "DICLAN" "EUPPER" "GALOBT" "LYSLAN"
## [129] "POTSIM" "SETPAR" "SOLMIS" "VIOSAG" "AGAFAS" "ANTNEG" "CROSAG" "LINSUL"
## [137] "RHUGLA" "FRAVIR" "GENPUB" "MUHGLA" "JUNBRA" "ULMSPP" "POLVER" "TRISPP"
## [145] "CORSPP" "POLSAN" "SILANT" "TRIPER" "RHUCOP"####### Remove Seed Rain Unknowns
# 1.) get a sum per species per transect of cover.
seed_rain_sum <- filter(seed_rain_sum, !grepl('UNK', SPP6))
seed_rain_sum <- filter(seed_rain_sum, !grepl('PANICUM?', SPP6))
seed_rain_sum <- filter(seed_rain_sum, !grepl('NONE', SPP6))
# Looks good
list(unique(seed_rain_sum$SPP6))## [[1]]
## [1] "ACAVIR" "AMBART" "ANAMIN" "BOLAST" "CERSPP" "CHAFAS" "CHEALB"
## [8] "CONCAN" "DIGISC" "ERISPP" "ERIVIL" "EUPSPP" "KUMSPP" "LEPVIR"
## [15] "MEDLUP" "OXADIL" "PENDIG" "PLAVIR" "PYCTEN" "RUDHIR" "SETSPP"
## [22] "SOLALT" "SORNUT" "SYMSPP" "THLARV" "TRIPER" "VEROSP" "ALOCAR"
## [29] "CORTRI" "CYPECH" "DICLAN" "ERYYUC" "HYPPUN" "MONFIS" "MYOVER"
## [36] "OENBIE" "RATPIN" "SCIGEO" "SOLRIG" "SPHOBT" "STRLEI" "TRIFLA"
## [43] "ACHMIL" "BARVUL" "FESARU" "GERCAR" "LESCUN" "LIAPYC" "PANCAP"
## [50] "SYMNOV" "BIDARI" "CARBUS" "DESSPP" "ERASPE" "GALAPA" "LESVIR"
## [57] "POAPRA" "RUDSUB" "VERHAS" "VERSPP" "ANDGER" "CARFES" "ECHCRU"
## [64] "EREHIE" "PLAOCC" "SILINT" "AGASPP" "BROJAP" "CAPBUR" "CIRALT"
## [71] "CORLAN" "MELSPP" "MOLVER" "SCHSCO" "SCLTRI" "CYPSTR" "FESPAR"
## [78] "JUNSPP" "PYCPIL" "TAROFF" "AMATUB" "LESCAP" "OENFIL" "KOEMAC"
## [85] "DAUCAR" "LEUVUL" "VULOCT" "CYPACU" "HORPUS" "SCIPEN" "ECHPAL"
## [92] "SOLNEM" "CORPAL" "LINSUL" "DIAARM" "HELMOL" "PERLON" "RUMCRI"
## [99] "BLECIL" "LYTALA" "SPOHET" "ELESPP1" "SPOCOM" "CARCEP" "AGRHYE"
## [106] "GENSPP" "HELHEL" "TRAOHI" "BAPALB" "CARANN" "IPOHED" "AMOCAN"
## [113] "SALAZU" "CARBIC" "EUPCOR" "EUTGYM" "GALOBT" "RUBSPP" "SETPAR"
## [120] "ELESPP" "LYSLAN" "POLSPP" "HYPHIR" "LYCAME" "SILANT" "CARDSP"
## [127] "CROSAG" "LOBSPI" "VIOSAG"####### Remove Seed Bank Unknowns
seed_bank_sum <- filter(seed_bank_sum, !grepl('UNK', SPP6))
# Looks good
list(sort(unique(seed_bank_sum$SPP6)))## [[1]]
## [1] "ABUTHE" "ACAVIR" "ACHMIL" "AGRHYE" "ALOCAR" "AMATUB" "AMBART" "ANDGER"
## [9] "BARVUL" "CARFES" "CARHIR" "CARPAR" "CERGLO" "CHAFAS" "CIRALT" "CONCAN"
## [17] "CORLAN" "CROSAG" "CYPACU" "CYPSPP" "CYPSTR" "DESSPP" "DICLAN" "DIGISC"
## [25] "DIGSAN" "DIPLAC" "ERASPE" "ERIANN" "ERISTR" "EUPHUM" "EUPSER" "EUTGYM"
## [33] "GERCAR" "HORPUS" "HYPPUN" "IPOHED" "JUNINT" "JUNMAR" "KUMSTI" "KUMSTR"
## [41] "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LYSLAN" "MEDLUP" "MELSPP" "MOLVER"
## [49] "OENSPP" "OXADIL" "PANCAP" "PANDIC" "PENDIG" "PERLON" "PLAOCC" "PLAPUS"
## [57] "PLAVIR" "POACHA" "POAPRA" "POPDEL" "POROLE" "PYCPIL" "PYCTEN" "RATPIN"
## [65] "RUDHIR" "RUMCRI" "SCHSCO" "SETFAB" "SETGLA" "SILANT" "SOLALT" "SOLCAR"
## [73] "SORNUT" "SPHOBT" "STRLEI" "SYMLAE" "SYMPIL" "TAROFF" "THLARV" "TOXRAD"
## [81] "TRAOHI" "TRIFLA" "TRIPER" "TRIREP" "VERHAS" "VERPER" "VERTHA"—–
Merge Spring and Fall cover
# Full join the two cover data sets
cover_merged <- full_join(spring_cover, fall_sum, by = c("Site", "Transect", "SPP6", "Scientific_Name", "Cover")) %>%
mutate_all(~replace(., is.na(.), 0))## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated# Make sure there is no weirdness in the join
#View(cover_merged)
### Group by and then select the SPP6 with the max cover between the two datasets
### Since Spring only had presence/absence we don't want to inflate the fall cover artificially if we saw that species again
cover_merged_max <- cover_merged %>%
group_by(Site, Transect, SPP6, Scientific_Name) %>%
summarise(Cover=max(Cover))## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.## Not exactly sure what to do about SYMSPP and KUMSPP (things we couldn't ID in the spring but could in the fall)
### Going to drop SYMSPP from the dataset since every plot that had SYMSPP in the spring had identified Symphyotrichum in the fall
cover_merged_max <- filter(cover_merged_max , !grepl('SYMSPP', SPP6))
### Checked KUMSPP only PFCA 2 - 3 and PFCA 2 - 9 did not see Kummerowia again in the fall; will drop KUMSPP for all other plots
### Removes all the PFCA 1 KUMSPPs
cover_merged_max <- cover_merged_max[!(cover_merged_max$SPP6%in%"KUMSPP" & cover_merged_max$Site%in% "PFCA 1"),]
### Removes all PFCA 2 KUMSPPs that were not in PFCA 2 - 3 and PFCA 2- 9
cover_merged_max <- cover_merged_max[!(cover_merged_max$SPP6%in%"KUMSPP" & cover_merged_max$Site%in% "PFCA 2" & cover_merged_max$Transect%in% c(2,4,5,8)),]
### Check all the species codes
sort(unique(cover_merged_max$SPP6))## [1] "ACAVIR" "ACERUB" "ACHMIL" "AGAFAS" "AGATEN" "AGRGIG" "AGRHYE" "ALOCAR"
## [9] "AMBART" "AMBPSI" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB"
## [17] "BAPAUS" "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD"
## [25] "CAPBUR" "CARBIC" "CARSPP" "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN"
## [33] "CORLAN" "CORPAL" "CORSPP" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR"
## [41] "DALPUR" "DESILL" "DESSES" "DESSPP" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL"
## [49] "ELETEN" "ELYCAN" "ELYVIR" "ERASPE" "ERISPP" "ERIVIL" "ERYYUC" "EUPALT"
## [57] "EUPCOR" "EUPPER" "EUPSER" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA"
## [65] "GALOBT" "GENAND" "GENPUB" "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP"
## [73] "HYPPUN" "JUNBRA" "JUNMAR" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "KUMSTI"
## [81] "KUMSTR" "LACSER" "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR"
## [89] "LIAPYC" "LIASPP" "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP"
## [97] "MELSPP" "MONFIS" "MUHGLA" "MUHSPP" "MYOVER" "OENBIE" "OENFIL" "OENSPP"
## [105] "OXADIL" "PANVIR" "PARINT" "PENDIG" "PLAMAJ" "PLAVIR" "POAPRA" "POLSAN"
## [113] "POLVER" "POTSIM" "PYCPIL" "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA"
## [121] "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO"
## [129] "SCISPP" "SCLTRI" "SENMAR" "SETFAB" "SETPAR" "SILANT" "SILINT" "SILLAC"
## [137] "SISSPP" "SOLALT" "SOLCAR" "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT"
## [145] "SPHOBT" "SPILAC" "SPOCOM" "SPOHET" "STRLEI" "SYMANO" "SYMERI" "SYMLAE"
## [153] "SYMLAN" "SYMLAT" "SYMNOV" "SYMOBL" "SYMOOL" "SYMORB" "SYMPIL" "SYMPRA"
## [161] "TAROFF" "THLARV" "TRAOHI" "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP"
## [169] "ULMPUM" "ULMSPP" "VERARV" "VERHEL" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT"
## [177] "ZIZAUR"—–
Merge Cover and Seed Rain
### Cover grouping so that it works with the seed rain dataset
# Code to lump certain species together
cover_to_merge_rain <- cover_merged_max
for(i in 1:nrow(cover_to_merge_rain)){
# Group all Symphyotrichum except SYMNOV
## SYMPIL
## SYMLAN
## SYMLAT
## SYMERI
## SYMANO
## SYMLAE
## SYMPRA
## SYMOBL
## SYMOOL
if(cover_to_merge_rain[i,3] == "SYMPIL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMLAN"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMLAT"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMERI"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMANO"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMLAE"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMPRA"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMOBL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
if(cover_to_merge_rain[i,3] == "SYMOOL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
# Group all Carex
## CARBIC
if(cover_to_merge_rain[i,3] == "CARBIC"){cover_to_merge_rain[i,3] <- "CARSPP"}
# Group all Kummerowia
## KUMSTI
## KUMSTR
if(cover_to_merge_rain[i,3] == "KUMSTI"){cover_to_merge_rain[i,3] <- "KUMSPP"}
if(cover_to_merge_rain[i,3] == "KUMSTR"){cover_to_merge_rain[i,3] <- "KUMSPP"}
# Group all Juncus
## JUNMAR
## JUNBBRA
if(cover_to_merge_rain[i,3] == "JUNMAR"){cover_to_merge_rain[i,3] <- "JUNSPP"}
if(cover_to_merge_rain[i,3] == "JUNBRA"){cover_to_merge_rain[i,3] <- "JUNSPP"}
# Group all Agalinis
## AGATEN
## AGAFAS
if(cover_to_merge_rain[i,3] == "AGATEN"){cover_to_merge_rain[i,3] <- "AGASPP"}
if(cover_to_merge_rain[i,3] == "AGAFAS"){cover_to_merge_rain[i,3] <- "AGASPP"}
# Group all Eleocharis
## ELETEN
if(cover_to_merge_rain[i,3] == "ELETEN"){cover_to_merge_rain[i,3] <- "ELESPP"}
# Group all Gentian
## GENPUB
## GENAND
if(cover_to_merge_rain[i,3] == "GENPUB"){cover_to_merge_rain[i,3] <- "GENSPP"}
if(cover_to_merge_rain[i,3] == "GENAND"){cover_to_merge_rain[i,3] <- "GENSPP"}
# Group all Setaria except Setaria parviflora
## SETGLA
## SETFAB
if(cover_to_merge_rain[i,3] == "SETGLA"){cover_to_merge_rain[i,3] <- "SETSPP"}
if(cover_to_merge_rain[i,3] == "SETFAB"){cover_to_merge_rain[i,3] <- "SETSPP"}
# Group all Veronica
## VERARV
if(cover_to_merge_rain[i,3] == "VERARV"){cover_to_merge_rain[i,3] <- "VEROSP"}
# Group all Desmodium
## DESSES
if(cover_to_merge_rain[i,3] == "DESSES"){cover_to_merge_rain[i,3] <- "DESSPP"}
# Group all Eupatorium
## EUPALT
## EUPSER
## EUPPER
if(cover_to_merge_rain[i,3] == "EUPALT"){cover_to_merge_rain[i,3] <- "EUPSPP"}
if(cover_to_merge_rain[i,3] == "EUPSER"){cover_to_merge_rain[i,3] <- "EUPSPP"}
if(cover_to_merge_rain[i,3] == "EUPPER"){cover_to_merge_rain[i,3] <- "EUPSPP"}
# Group all Polygala
## POLVER
## POLSAN
if(cover_to_merge_rain[i,3] == "POLVER"){cover_to_merge_rain[i,3] <- "POLSPP"}
if(cover_to_merge_rain[i,3] == "POLSAN"){cover_to_merge_rain[i,3] <- "POLSPP"}
}
### Check that it changed the codes
sort(unique(cover_to_merge_rain$SPP6))## [1] "ACAVIR" "ACERUB" "ACHMIL" "AGASPP" "AGRGIG" "AGRHYE" "ALOCAR" "AMBART"
## [9] "AMBPSI" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB" "BAPAUS"
## [17] "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD" "CAPBUR"
## [25] "CARSPP" "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN" "CORLAN" "CORPAL"
## [33] "CORSPP" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR" "DALPUR" "DESILL"
## [41] "DESSPP" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL" "ELESPP" "ELYCAN" "ELYVIR"
## [49] "ERASPE" "ERISPP" "ERIVIL" "ERYYUC" "EUPCOR" "EUPSPP" "EUTGYM" "FESARU"
## [57] "FESPAR" "FRAVIR" "GALAPA" "GALOBT" "GENSPP" "GERCAR" "HELFLE" "HELHEL"
## [65] "HELMOL" "HELSPP" "HYPPUN" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER"
## [73] "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR" "LIAPYC" "LIASPP"
## [81] "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP" "MELSPP" "MONFIS"
## [89] "MUHGLA" "MUHSPP" "MYOVER" "OENBIE" "OENFIL" "OENSPP" "OXADIL" "PANVIR"
## [97] "PARINT" "PENDIG" "PLAMAJ" "PLAVIR" "POAPRA" "POLSPP" "POTSIM" "PYCPIL"
## [105] "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA" "ROSCAR" "ROSSPP" "RUBSPP"
## [113] "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO" "SCISPP" "SCLTRI" "SENMAR"
## [121] "SETPAR" "SETSPP" "SILANT" "SILINT" "SILLAC" "SISSPP" "SOLALT" "SOLCAR"
## [129] "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT" "SPHOBT" "SPILAC" "SPOCOM"
## [137] "SPOHET" "STRLEI" "SYMNOV" "SYMORB" "SYMSPP" "TAROFF" "THLARV" "TRAOHI"
## [145] "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP" "ULMPUM" "ULMSPP" "VERHEL"
## [153] "VEROSP" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT" "ZIZAUR"### Sum everything by species by transect again
cover_to_merge_rain <- cover_to_merge_rain %>%
group_by(Site, Transect, SPP6) %>%
summarize(Cover = sum(Cover))## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.### Seed grouping so that it works with the cover dataset
seed_rain_names <- left_join(seed_rain_sum, traits)## Joining with `by = join_by(SPP6)`rain_to_merge_cover <- seed_rain_names[,c(1,2,3,4,5)]
for(i in 1:nrow(rain_to_merge_cover)){
# Group all Carex
## CARBUS
## CARBIC
## CARANN
## CARFES
## CARCEP
if(rain_to_merge_cover[i,3] == "CARBUS"){rain_to_merge_cover[i,3] <- "CARSPP"}
if(rain_to_merge_cover[i,3] == "CARBIC"){rain_to_merge_cover[i,3] <- "CARSPP"}
if(rain_to_merge_cover[i,3] == "CARANN"){rain_to_merge_cover[i,3] <- "CARSPP"}
if(rain_to_merge_cover[i,3] == "CARFES"){rain_to_merge_cover[i,3] <- "CARSPP"}
if(rain_to_merge_cover[i,3] == "CARCEP"){rain_to_merge_cover[i,3] <- "CARSPP"}
# Group all Eleocharis
## ELESPP1
if(rain_to_merge_cover[i,3] == "ELESPP1"){rain_to_merge_cover[i,3] <- "ELESPP"}
# Group all Scirpus
if(rain_to_merge_cover[i,3] == "SCIGEO"){rain_to_merge_cover[i,3] <- "SCISPP"}
if(rain_to_merge_cover[i,3] == "SCIPEN"){rain_to_merge_cover[i,3] <- "SCISPP"}
}
### Sum everything by species by transect again
rain_to_merge_cover <- rain_to_merge_cover %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seeds = sum(Number_Seeds))## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.# Looks good
sort(unique(rain_to_merge_cover$SPP6))## [1] "ACAVIR" "ACHMIL" "AGASPP" "AGRHYE" "ALOCAR" "AMATUB" "AMBART" "AMOCAN"
## [9] "ANAMIN" "ANDGER" "BAPALB" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP"
## [17] "CAPBUR" "CARDSP" "CARSPP" "CERSPP" "CHAFAS" "CHEALB" "CIRALT" "CONCAN"
## [25] "CORLAN" "CORPAL" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR" "DAUCAR"
## [33] "DESSPP" "DIAARM" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL" "ELESPP" "ERASPE"
## [41] "EREHIE" "ERISPP" "ERIVIL" "ERYYUC" "EUPCOR" "EUPSPP" "EUTGYM" "FESARU"
## [49] "FESPAR" "GALAPA" "GALOBT" "GENSPP" "GERCAR" "HELHEL" "HELMOL" "HORPUS"
## [57] "HYPHIR" "HYPPUN" "IPOHED" "JUNSPP" "KOEMAC" "KUMSPP" "LEPVIR" "LESCAP"
## [65] "LESCUN" "LESVIR" "LEUVUL" "LIAPYC" "LINSUL" "LOBSPI" "LYCAME" "LYSLAN"
## [73] "LYTALA" "MEDLUP" "MELSPP" "MOLVER" "MONFIS" "MYOVER" "OENBIE" "OENFIL"
## [81] "OXADIL" "PANCAP" "PENDIG" "PERLON" "PLAOCC" "PLAVIR" "POAPRA" "POLSPP"
## [89] "PYCPIL" "PYCTEN" "RATPIN" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU"
## [97] "SCHSCO" "SCISPP" "SCLTRI" "SETPAR" "SETSPP" "SILANT" "SILINT" "SOLALT"
## [105] "SOLNEM" "SOLRIG" "SORNUT" "SPHOBT" "SPOCOM" "SPOHET" "STRLEI" "SYMNOV"
## [113] "SYMSPP" "TAROFF" "THLARV" "TRAOHI" "TRIFLA" "TRIPER" "VERHAS" "VEROSP"
## [121] "VERSPP" "VIOSAG" "VULOCT"# Drop the scientific names they are causing trouble with the join
cover_to_merge_rain <- cover_to_merge_rain
cover_rain <- full_join(cover_to_merge_rain, rain_to_merge_cover) %>%
mutate_all(~replace(., is.na(.), 0))## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:# Drop TP transect 7 because of missing data
cover_rain <- cover_rain[!(cover_rain$Site%in%"TP" & cover_rain$Transect%in% "7"),]—–
Merge Cover and Seed Bank
### Cover grouping so that it works with the seed bank dataset
# Code to lump certain species together
cover_to_merge_bank <- cover_merged_max
for(i in 1:nrow(cover_to_merge_bank)){
# Group all Desmodium together
## DESSES
if(cover_to_merge_bank[i,3] == "DESSES"){cover_to_merge_bank[i,3] <- "DESSPP"}
# Group all Cyperus together
## CYPSTR
## CYPECH
if(cover_to_merge_bank[i,3] == "CYPSTR"){cover_to_merge_bank[i,3] <- "CYPSPP"}
if(cover_to_merge_bank[i,3] == "CYPECH"){cover_to_merge_bank[i,3] <- "CYPSPP"}
# Group all Oenothera together
## OENFIL
## OENBIE
if(cover_to_merge_bank[i,3] == "OENFIL"){cover_to_merge_bank[i,3] <- "OENSPP"}
if(cover_to_merge_bank[i,3] == "OENBIE"){cover_to_merge_bank[i,3] <- "OENSPP"}
}
cover_to_merge_bank <- cover_to_merge_bank %>%
group_by(Site, Transect, SPP6) %>%
summarize(Cover = sum(Cover))## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.### Seed Bank grouping so that it works with the cover dataset
bank_to_merge_cover <- seed_bank_sum
# Code to lump certain species together
for(i in 1:nrow(bank_to_merge_cover)){
# Group all Cerastium together
## CERGLO
if(bank_to_merge_cover[i,3] == "CERGLO"){bank_to_merge_cover[i,3] <- "CERSPP"}
# Group all Cyperus together
## DESSES
if(bank_to_merge_cover[i,3] == "CYPSTR"){bank_to_merge_cover[i,3] <- "CYPSPP"}
# Group all the Carex together
if(bank_to_merge_cover[i,3] == "CARFES"){bank_to_merge_cover[i,3] <- "CARSPP"}
# Group all Erigeron together
## ERIANN
## ERISTR
if(bank_to_merge_cover[i,3] == "ERIANN"){bank_to_merge_cover[i,3] <- "ERISPP"}
if(bank_to_merge_cover[i,3] == "ERISTR"){bank_to_merge_cover[i,3] <- "ERISPP"}
# Change Juncus interior to Juncus sp. (since we were able to ID JUNMAR in the cover but not JUNINT)
if(bank_to_merge_cover[i,3] == "JUNINT"){bank_to_merge_cover[i,3] <- "JUNSPP"}
}
bank_to_merge_cover <- bank_to_merge_cover %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seedlings = sum(Number_Seedlings))## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.# Drop the scientific names they are causing trouble with the join
cover_bank <- full_join(cover_to_merge_bank, bank_to_merge_cover) %>%
mutate_all(~replace(., is.na(.), 0))## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:# Drop TP transect 7 because of missing fall cover data
cover_bank <- cover_bank[!(cover_bank$Site%in%"TP" & cover_bank$Transect%in% "7"),]—–
Merge Seed Rain and Seed Bank
### Seed Bank grouping so that it works with the cover dataset
bank_to_merge_rain <- seed_bank_sum
# Code to lump certain species together
for(i in 1:nrow(bank_to_merge_rain)){
# Group all Cerastium together
## CERGLO
if(bank_to_merge_rain[i,3] == "CERGLO"){bank_to_merge_rain[i,3] <- "CERSPP"}
# Group all Cyperus together
## CYPSTR
if(bank_to_merge_rain[i,3] == "CYPSTR"){bank_to_merge_rain[i,3] <- "CYPSPP"}
# Group all Erigeron together
## ERIANN
## ERISTR
if(bank_to_merge_rain[i,3] == "ERIANN"){bank_to_merge_rain[i,3] <- "ERISPP"}
if(bank_to_merge_rain[i,3] == "ERISTR"){bank_to_merge_rain[i,3] <- "ERISPP"}
# Group all the Juncus together
## JUNINT
## JUNMAR
if(bank_to_merge_rain[i,3] == "JUNINT"){bank_to_merge_rain[i,3] <- "JUNSPP"}
if(bank_to_merge_rain[i,3] == "JUNMAR"){bank_to_merge_rain[i,3] <- "JUNSPP"}
# Group all the Setaria together
if(bank_to_merge_rain[i,3] == "SETFAB"){bank_to_merge_rain[i,3] <- "SETSPP"}
if(bank_to_merge_rain[i,3] == "SETGLA"){bank_to_merge_rain[i,3] <- "SETSPP"}
# Group all the Symphyotrichum together
## SYMPIL
## SYMLAE
if(bank_to_merge_rain[i,3] == "SYMPIL"){bank_to_merge_rain[i,3] <- "SYMSPP"}
if(bank_to_merge_rain[i,3] == "SYMLAE"){bank_to_merge_rain[i,3] <- "SYMSPP"}
if(bank_to_merge_rain[i,3] == "SIMPIL"){bank_to_merge_rain[i,3] <- "SYMSPP"}
# Group all the Kummerowia together
## KUMSTR
## KUMSTI
if(bank_to_merge_rain[i,3] == "KUMSTR"){bank_to_merge_rain[i,3] <- "KUMSPP"}
if(bank_to_merge_rain[i,3] == "KUMSTI"){bank_to_merge_rain[i,3] <- "KUMSPP"}
# Group all the Veronica together
## VERPER
if(bank_to_merge_rain[i,3] == "VERPER"){bank_to_merge_rain[i,3] <- "VEROSP"}
# Group all the Eupatorium together
## EUPSER
if(bank_to_merge_rain[i,3] == "EUPSER"){bank_to_merge_rain[i,3] <- "EUPSPP"}
# Merge all the Cardamine together
## CARHIR
## CARPAR
if(bank_to_merge_rain[i,3] == "CARHIR"){bank_to_merge_rain[i,3] <- "CARDSP"}
if(bank_to_merge_rain[i,3] == "CARPAR"){bank_to_merge_rain[i,3] <- "CARDSP"}
}
bank_to_merge_rain <- bank_to_merge_rain %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seedlings = sum(Number_Seedlings))## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.rain_to_merge_bank <- seed_rain_sum
for(i in 1:nrow(rain_to_merge_bank)){
# Group all Cyperus together
## CYPSTR
## CYPECH
if(rain_to_merge_bank[i,3] == "CYPSTR"){rain_to_merge_bank[i,3] <- "CYPSPP"}
if(rain_to_merge_bank[i,3] == "CYPECH"){rain_to_merge_bank[i,3] <- "CYPSPP"}
# Group all Oenothera together
## OENFIL
## OENBIE
if(rain_to_merge_bank[i,3] == "OENFIL"){rain_to_merge_bank[i,3] <- "OENSPP"}
if(rain_to_merge_bank[i,3] == "OENBIE"){rain_to_merge_bank[i,3] <- "OENSPP"}
}
rain_to_merge_bank <- rain_to_merge_bank %>%
group_by(Site, Transect, SPP6) %>%
summarize(Number_Seeds = sum(Number_Seeds))## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.rain_bank <- full_join(bank_to_merge_rain, rain_to_merge_bank) %>%
mutate_all(~replace(., is.na(.), 0))## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:—-
Seed Rain and Cover Exploratory Analysis
cover_rain_sum <- cover_rain %>%
group_by(SPP6) %>%
summarize(tot.cover = sum(Cover), tot.seeds = sum(Number_Seeds))Wasn’t in the cover at all but we caught seeds
- Amaranthus tuberculatus
- Cardamine spp.
- Chenopodium album
- Daucus carota
- Dianthus armeria
- Erechtites hieracifolia
- Hordeum pusillum
- Ipomoea hederacea
- Leucanthemum vulgare
- Mollugo verticillata
- Panicum capillare
- Persicaria longiseta
- Platanus occidentalis
- Verbena hastata
Didn’t catch seeds but was in the cover
Notable (> 10% cover observed)
- Baptisia australis
- Potentilla simplex
- Silphium laciniatum
- Rosa carolina
- Fragraria virginica
- Campsis radicans
- Parthenium integrifolium
- Rhus copallinum
- Solanum carolinense
- Helianthus spp (the pressed one)
- Baptisia bracteata
- Antennaria neglecta
Minor ( 1% < cover < 10% observed)
- Elymus virginica
- Rhus glabra
- Trifolium campestre
- Verbesina helianthoides
- Solidago speciosa
- Lespedeza procombens
- Liatris spp.
- Ranunculus abortivus
- Juniperus virginiana
- Trifolium spp.
- Ulmus spp.
- Agrostis gigantea
- Lactuca serriola
- Symphoricarpos orbiculatus
- Ulmus pumila
- Acer rubrum
- Comandra umbellata
- Dalea purpurea
- Muhlenbergia spp.
- Plantago major
- Rosa spp.
- Spiranthes lacera
- Trifolium repens
- Viola spp.
- Zizea aurea
Trace (cover 1% observed)
- Ambrosia psilostachya
- Asclepias syriaca
- Desmanthus illinoensis
- Elymus canadense
- Helenium flexuosum
- Lepidium densiflorum
- Liatris squarrosa
- Muhlenbergia glabrifolia
- Oenothera spp. (not biennis)
- Panicum virgatum
- Senna marilandica
- Sisyrinchium spp.
- Solidago missouriensis.
Dissimilarity analysis
## # A tibble: 10,959 × 5
## # Groups: Site, Transect [39]
## Site Transect SPP6 Abundance Type
## <chr> <fct> <chr> <dbl> <chr>
## 1 PFCA 1 1 ACAVIR 11.5 Aboveground
## 2 PFCA 1 10 ACAVIR 4 Aboveground
## 3 PFCA 1 2 ACAVIR 4.5 Aboveground
## 4 PFCA 1 3 ACAVIR 2 Aboveground
## 5 PFCA 1 4 ACAVIR 4 Aboveground
## 6 PFCA 1 5 ACAVIR 3 Aboveground
## 7 PFCA 1 6 ACAVIR 2.5 Aboveground
## 8 PFCA 1 7 ACAVIR 3 Aboveground
## 9 PFCA 1 8 ACAVIR 1 Aboveground
## 10 PFCA 1 9 ACAVIR 5 Aboveground
## # ℹ 10,949 more rows# Code obtained and modified from Lauren and Anu's seed rain vs. vegetation paper
##### run dissimilarity loops
Cover_rain_dis_results <-matrix(nrow=0,ncol=7)
sites <-as.vector(unique(cover_rain_dis$Site))
for(s in 1:length(sites)){
temp <-subset(cover_rain_dis, Site==sites[s])
temp$SPP6 <- factor(as.character(temp$SPP6))
transects <-as.vector(unique(temp$Transect))
#for(b in 1:length(blocks)){
#temp_b<-subset(temp, block==blocks[b])
#plots<-as.vector(unique(temp_b$plot))
for(t in 1:length(transects)){
temp_t <- subset(temp, Transect==transects[t])
#wide form - with a row for above and belowground
plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
#relative abundances
trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"tot"))
names(trans_tot)[1]<-"trt"
trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
names(trans_pa)[1]<-"trt"
trans_int <- cbind(plot_cast$Type,trunc(plot_cast[,-1])) #integers only
names(trans_int)[1]<-"trt"
trans_int_norm <- cbind(plot_cast$Type,trunc(decostand(plot_cast[,-1], "normalize")*100)) #integers only and normalize
names(trans_int)[1]<-"trt"
distance_bray<-vegdist(trans_tot[,-1], method = "bray")
distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")
distance_kulczynski<-vegdist(trans_tot[,-1], method = "kulczynski")
distance_cao<-vegdist(trans_int_norm[,-1], method = "cao")
distance_morisita<-vegdist(trans_int[,-1], method = "morisita")
new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1],
distance_cao[1], distance_morisita[1], distance_kulczynski[1])
Cover_rain_dis_results <-rbind(Cover_rain_dis_results, new.row)
}
print(s/length(sites))
}## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1row.names(Cover_rain_dis_results)<-NULL
Cover_rain_dis_results<- data.frame(Cover_rain_dis_results)
names(Cover_rain_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard",
"dissimilarity_cao", "dissimilarity_morisita", "dissimilarity_kulczynski")##
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = Cover_rain_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.267522 -0.063391 -0.002319 0.088438 0.190470
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.680648 0.038164 17.835 < 2e-16 ***
## SitePFCA 1 -0.211837 0.052605 -4.027 0.000289 ***
## SitePFCA 2 -0.049672 0.052605 -0.944 0.351519
## SitePFCA 3 0.006873 0.052605 0.131 0.896792
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1145 on 35 degrees of freedom
## Multiple R-squared: 0.4014, Adjusted R-squared: 0.3501
## F-statistic: 7.823 on 3 and 35 DF, p-value: 0.0003994## Anova Table (Type III tests)
##
## Response: dissimilarity_bray
## Sum Sq Df F value Pr(>F)
## (Intercept) 4.1695 1 318.0843 < 2.2e-16 ***
## Site 0.3076 3 7.8231 0.0003994 ***
## Residuals 0.4588 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Cover_rain.mod.bray)
##
## $Site
## diff lwr upr p adj
## PFCA 1-TP -0.211837483 -0.35370841 -0.06996655 0.0015759
## PFCA 2-TP -0.049671993 -0.19154292 0.09219894 0.7813451
## PFCA 3-TP 0.006873424 -0.13499751 0.14874435 0.9991874
## PFCA 2-PFCA 1 0.162165490 0.02407847 0.30025251 0.0160642
## PFCA 3-PFCA 1 0.218710907 0.08062388 0.35679793 0.0007822
## PFCA 3-PFCA 2 0.056545416 -0.08154161 0.19463244 0.6892609##
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = Cover_rain_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.104048 -0.037405 -0.006949 0.037294 0.192794
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.54155 0.02038 26.578 < 2e-16 ***
## SitePFCA 1 -0.06193 0.02809 -2.205 0.03412 *
## SitePFCA 2 -0.07926 0.02809 -2.822 0.00782 **
## SitePFCA 3 0.04722 0.02809 1.681 0.10162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06113 on 35 degrees of freedom
## Multiple R-squared: 0.4353, Adjusted R-squared: 0.3869
## F-statistic: 8.994 on 3 and 35 DF, p-value: 0.0001491## Anova Table (Type III tests)
##
## Response: dissimilarity_jaccard
## Sum Sq Df F value Pr(>F)
## (Intercept) 2.63947 1 706.3882 < 2.2e-16 ***
## Site 0.10082 3 8.9943 0.0001491 ***
## Residuals 0.13078 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Cover_rain.mod.jac)
##
## $Site
## diff lwr upr p adj
## PFCA 1-TP -0.06192816 -0.13767379 0.013817463 0.1417750
## PFCA 2-TP -0.07925589 -0.15500152 -0.003510265 0.0374144
## PFCA 3-TP 0.04721930 -0.02852633 0.122964926 0.3485548
## PFCA 2-PFCA 1 -0.01732773 -0.09105311 0.056397651 0.9203909
## PFCA 3-PFCA 1 0.10914746 0.03542208 0.182872842 0.0017366
## PFCA 3-PFCA 2 0.12647519 0.05274981 0.200200570 0.0002772Fig. - Bray-Curtis dissimilarity
## 95% confidence intervals
emmeans::emmeans(Cover_rain.mod.bray, "Site", type = "response")## Site emmean SE df lower.CL upper.CL
## TP 0.681 0.0382 35 0.603 0.758
## PFCA 1 0.469 0.0362 35 0.395 0.542
## PFCA 2 0.631 0.0362 35 0.557 0.704
## PFCA 3 0.688 0.0362 35 0.614 0.761
##
## Confidence level used: 0.95## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.469, 0.631, 0.688, 0.681)
lower <- c(0.395, 0.557, 0.614, 0.603)
upper <- c(0.542, 0.704, 0.761, 0.758)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Cover_rain_dissimilarity <- data.frame(site, response, lower, upper)
Fig. - Jaccard dissimilarity figure
## 95% confidence intervals
emmeans::emmeans(Cover_rain.mod.jac, "Site", type = "response")## Site emmean SE df lower.CL upper.CL
## TP 0.542 0.0204 35 0.500 0.583
## PFCA 1 0.480 0.0193 35 0.440 0.519
## PFCA 2 0.462 0.0193 35 0.423 0.502
## PFCA 3 0.589 0.0193 35 0.550 0.628
##
## Confidence level used: 0.95## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.480, 0.462, 0.589, 0.542)
lower <- c(0.440, .423, 0.550, 0.500)
upper <- c(0.519, 0.502, 0.628, 0.583)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Cover_rain_jac_dissimilarity <- data.frame(site, response, lower, upper)
Immigration
Proportion of new species in the seed rain not seen in local cover
cover_rain_dis_wide <- cover_rain_dis %>%
spread(SPP6, Abundance) %>%
mutate_all(~replace(., is.na(.), 0)) ## `mutate_all()` ignored the following grouping variables:
## • Columns `Site`, `Transect`
## ℹ Use `mutate_at(df, vars(-group_cols()), myoperation)` to silence the message.cover_rain_dis_long <- cover_rain_dis_wide %>%
gather(key = "SPP6", value = "Abundance", 4:177) %>%
spread(key = Type, Abundance)
# Get rid of columns that both have zeros
cover_rain_dis_long_unique <- cover_rain_dis_long[!(cover_rain_dis_long$Aboveground == 0 & cover_rain_dis_long$Seed_Rain == 0),]
# 0 = no seeds/seedlings and 1 = seeds/seedlings
cover_rain_dis_long_unique <- cover_rain_dis_long_unique %>%
mutate(Aboveground = if_else(Aboveground > 0, 1, 0)) %>%
mutate(Seed_Rain = if_else(Seed_Rain > 0, 1, 0))
# Make a new column that adds both the seed rain and aboveground presence/absence columns. Then remove rows that = 2 (shares both species)
cover_rain_dis_long_unique$shared <- cover_rain_dis_long_unique$Aboveground + cover_rain_dis_long_unique$Seed_Rain
cover_rain_dis_long_unique2 <- cover_rain_dis_long_unique[!(cover_rain_dis_long_unique$shared == 2),]
cover_rain_dis_long_unique2 <- cover_rain_dis_long_unique2[, -6]
# Count unique taxa per site and transect for the seed bank compared to the seed rain
Unique_species_aboveground_com_rain <- cover_rain_dis_long_unique2[, -5]
Unique_species_aboveground_com_rain <- Unique_species_aboveground_com_rain %>%
filter(Aboveground != 0 ) %>%
group_by(Site, Transect) %>%
summarize(New_Aboveground = length(unique(SPP6)))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.# Count unique taxa per site and transect for the seed rain compared to the seed bank
Unique_species_rain_com_aboveground <- cover_rain_dis_long_unique2[, -4]
Unique_species_rain_com_aboveground <- Unique_species_rain_com_aboveground %>%
filter(Seed_Rain != 0 ) %>%
group_by(Site, Transect) %>%
summarize(New_Rain = length(unique(SPP6)))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.# Count the number of shared taxa between communities
cover_rain_dis_long_shared <- cover_rain_dis_long_unique %>%
filter(shared == 2) %>%
group_by(Site, Transect) %>%
summarize(Shared = length(unique(SPP6)))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.## Merge them together
Species_counts_cover_rain <- full_join(Unique_species_aboveground_com_rain, Unique_species_rain_com_aboveground) %>%
mutate_all(~replace(., is.na(.), 0)) ## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.Species_counts_cover_rain <- full_join(Species_counts_cover_rain, cover_rain_dis_long_shared) %>%
mutate_all(~replace(., is.na(.), 0)) ## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.Species_counts_cover_rain$Tot_Richness <- Species_counts_cover_rain$New_Aboveground + Species_counts_cover_rain$New_Rain + Species_counts_cover_rain$Shared
# Count the number of species in the seed rain for each transect
Species_counts_cover_rain$Tot_Rain <- Species_counts_cover_rain$New_Rain + Species_counts_cover_rain$Shared
# Count the number of species in the seed bank for each transect
Species_counts_cover_rain$Tot_Aboveground <- Species_counts_cover_rain$New_Aboveground + Species_counts_cover_rain$Shared
## Immigrants
Species_counts_cover_rain$Prop_Immigrant <- Species_counts_cover_rain$New_Rain / Species_counts_cover_rain$Tot_Abovegroundfavstats(Species_counts_cover_rain$Prop_Immigrant~Species_counts_cover_rain$Site)## Species_counts_cover_rain$Site min Q1 median Q3
## 1 PFCA 1 0.2325581 0.2454268 0.3009259 0.4000000
## 2 PFCA 2 0.1621622 0.2231378 0.2889714 0.3224343
## 3 PFCA 3 0.1081081 0.2181818 0.3042017 0.3495879
## 4 TP 0.1851852 0.2162162 0.2592593 0.3103448
## max mean sd n missing
## 1 0.6190476 0.3420640 0.12881067 10 0
## 2 0.3555556 0.2738325 0.06337656 10 0
## 3 0.4736842 0.2832612 0.11509477 10 0
## 4 0.5161290 0.2849422 0.10421660 9 0# No difference in proportion of immigrant species
## *Note* better modeled as a binomial regression?
Species_counts_cover_rain.mod <- lm(Prop_Immigrant~Site, data =Species_counts_cover_rain)
summary(Species_counts_cover_rain.mod)##
## Call:
## lm(formula = Prop_Immigrant ~ Site, data = Species_counts_cover_rain)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.17515 -0.07599 -0.01296 0.05072 0.27698
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.34206 0.03345 10.226 4.7e-12 ***
## SitePFCA 2 -0.06823 0.04730 -1.442 0.158
## SitePFCA 3 -0.05880 0.04730 -1.243 0.222
## SiteTP -0.05712 0.04860 -1.175 0.248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1058 on 35 degrees of freedom
## Multiple R-squared: 0.06862, Adjusted R-squared: -0.01121
## F-statistic: 0.8596 on 3 and 35 DF, p-value: 0.4711Anova(Species_counts_cover_rain.mod, type = "III")## Anova Table (Type III tests)
##
## Response: Prop_Immigrant
## Sum Sq Df F value Pr(>F)
## (Intercept) 1.17008 1 104.5809 4.703e-12 ***
## Site 0.02885 3 0.8596 0.4711
## Residuals 0.39159 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1prop_new_species_cover_rain.plot <- ggplot()+
geom_jitter(data=Species_counts_cover_rain, aes(x = Site, y = Prop_Immigrant, color = Site), show.legend = FALSE, width = .2, size = 2.5) +
theme_classic()+
labs(x = "", y = "Prop. of new species in seed rain")+
scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"))+
ylim(-0.05,.75)
prop_new_species_cover_rain.plot
Proportion of seeds in the seed rain not seen in local cover
names(cover_rain_dis_long) <- c("Site", "Transect", "SPP6", "Aboveground_Abundance", "Seed_Rain_Abundance")
cover_rain_dis_long_unique_abundance <- full_join(cover_rain_dis_long_unique, cover_rain_dis_long)## Joining with `by = join_by(Site, Transect, SPP6)`cover_rain_dis_long_shared_abundance <- cover_rain_dis_long_unique_abundance %>%
filter(shared == 2) %>%
group_by(Site, Transect) %>%
summarize(Shared_Seeds = sum(Seed_Rain_Abundance))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.cover_rain_dis_long_immigrant_abundance <- cover_rain_dis_long_unique_abundance %>%
filter(shared == 1) %>%
group_by(Site, Transect) %>%
summarize(Immigrant_Seeds = sum(Seed_Rain_Abundance))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.Prop_Immigrant_Seeds <- full_join(cover_rain_dis_long_shared_abundance, cover_rain_dis_long_immigrant_abundance)## Joining with `by = join_by(Site, Transect)`Prop_Immigrant_Seeds$Total_Seeds <- Prop_Immigrant_Seeds$Shared_Seeds + Prop_Immigrant_Seeds$Immigrant_Seeds
Prop_Immigrant_Seeds$Prop <- Prop_Immigrant_Seeds$Immigrant_Seeds / Prop_Immigrant_Seeds$Total_Seedsprop_immigrant_seeds_cover_rain.plot <- ggplot()+
geom_jitter(data=Prop_Immigrant_Seeds, aes(x = Site, y = Prop, color = Site), show.legend = FALSE, width = .2, size = 2.5) +
theme_classic()+
labs(x = "", y = "Prop. of new seeds in seed rain", title = "Aboveground - Seed Rain")+
scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), plot.title = element_text(hjust = 0.5))+
ylim(-0.05,.5)
# The two "outlier" points are because we caught a lot of Solidago altissima and Barbarea vulgaris at these points but they weren't in the immediate cover.
prop_immigrant_seeds_cover_rain.plot
immigrant_resident_mass <- cover_rain_dis_long_unique %>%
filter(Seed_Rain == 1 & Aboveground == 0)
immigrant_resident_mass_shared <- cover_rain_dis_long_unique %>%
filter(Seed_Rain == 1 & Aboveground == 1)
immigrant_resident_mass <- full_join(immigrant_resident_mass, immigrant_resident_mass_shared)## Joining with `by = join_by(Site, Transect, SPP6, Aboveground, Seed_Rain,
## shared)`immigrant_resident_mass <- left_join(immigrant_resident_mass , traits)## Joining with `by = join_by(SPP6)`ggplot(aes(x = as.factor(shared), y = Mean_1_Seed_Mass_g*1000), data = immigrant_resident_mass)+
geom_boxplot() +
scale_y_continuous(trans='log10') +
theme_classic()+
labs(x = "", y = "Seed mass (mg) [log 10]")+
scale_x_discrete(labels = c("Immigrant", "Resident"), breaks = c("1", "2"))+
theme(text=element_text(size=18), legend.key.size=unit(0.25, "cm"))## Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
test_mod <- lm(log10(Mean_1_Seed_Mass_g*1000) ~ as.factor(shared), data = immigrant_resident_mass)
anova(test_mod)## Analysis of Variance Table
##
## Response: log10(Mean_1_Seed_Mass_g * 1000)
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(shared) 1 6.34 6.336 14.176 0.0001741 ***
## Residuals 1279 571.66 0.447
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ggplot(aes(x = as.factor(shared), y = Mean_1_Seed_Mass_g*1000), data = immigrant_resident_mass)+
geom_boxplot() +
scale_y_continuous(trans='log10') +
theme_classic()+
labs(x = "", y = "Seed mass (mg) [log 10]")+
scale_x_discrete(labels = c("Immigrant", "Resident"), breaks = c("1", "2"))+
theme(text=element_text(size=18), legend.key.size=unit(0.25, "cm"))## Warning: Removed 2 rows containing non-finite values (`stat_boxplot()`).
Cover vs. Seed Rain - Species-level and traits
By Site
#cover_rain$Site <- factor(cover_rain$Site, levels = c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))
site.labs <- c("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant")
ggplot(cover_rain, aes(x = Cover+1, y = Number_Seeds+1))+
geom_point()+
scale_x_continuous(trans='log10') +
scale_y_continuous(trans='log10') +
facet_grid(cols = vars(Site), labeller = labeller(Site = site.labs))+
geom_smooth(method = "lm") +
labs(x = "Veg. Cover [Log 10 + 1])", y = "No. Seeds [Log 10 + 1]")+
theme(axis.text = element_text( size = 12),
axis.text.x = element_text( size = 12),
axis.title = element_text( size = 16),
legend.position = "none",
strip.text = element_text(size = 16))## `geom_smooth()` using formula = 'y ~ x'
Residual plots
# # Remove everything that we didn't see both cover and seeds for at a transect
# ## Removing the "rare" occurrences (cover has to be greater than 1 and number of seeds greater than 10) removes the streaking apparent in the residuals
#
# cover_rain_filtered <- cover_rain %>%
# filter(Cover > 0) %>%
# filter(Number_Seeds > 0)# ggplot(cover_rain_filtered, aes(x = Cover, y = Number_Seeds, color = Site))+
# geom_point()+
# scale_x_continuous(trans='log10') +
# scale_y_continuous(trans='log10') +
# geom_smooth(method = "lm", alpha = 0.2) +
# labs(y = "No. Seeds [Log 10 + 1]", x = "Veg. Cover [Log 10 + 1]", title = "Seed rain vs cover")+
# scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
# theme_classic() +
# theme(axis.text = element_text( size = 12),
# axis.text.x = element_text( size = 12),
# axis.title = element_text( size = 16),
# legend.position = "none",
# strip.text = element_text(size = 16))Fitting the model
Note I tried other distributional families (gamma, inverse guassian, tweedie, negative binomial) and logging both variables and using a Normal distribution seems to be the best option when looking at the residual vs. fits plot. You do get weird streaks in the residuals though because there were so many instances of us seeing 1 or 2 seeds/cover.
#
# # Fit a nb mixed-model
#
# mod1 <- glmer.nb(Number_Seeds~log10(Cover+1)+log10(Cover+1):Site+(1|Site), data = cover_rain_filtered)
# summary(mod1)
# r.squaredGLMM(mod1)
#
#
# check_model(mod1)
#
# plot(mod1)
#
# check_overdispersion(mod1)
# Histograms of the predictor and response variables after logging show why we get the streaks. They still aren’t particularly Normal after logging.
# p <-ggplot( cover_rain_filtered, aes(x=(log10(prop.seeds+0.01)))) +
# geom_histogram(color="black", fill="white")+
# labs(x = "Log(no. seeds)", y = "Frequency")
# p
#
#
# q <-ggplot(cover_rain_filtered, aes(x=log10(prop.cover+0.01))) +
# geom_histogram(color="black", fill="white")+
# labs(x = "Log(total cover)", y = "Frequency")
# q Plotting the species and their associated residuals shows which species were over and under represented in the seed rain based on local vegetative cover. It does not show instances where we caught immigrating seeds (no local cover) or species with unsuccessful reproduction (no captured seeds).
# resid.cover.rain <- residuals(mod1)
# SPP6 <- cover_rain_filtered$SPP6
# Site <- cover_rain_filtered$Site
#
# resid.species.all.cover.rain <- data.frame(SPP6, resid.cover.rain, Site)
# resid.specieslall.cover.rain <- resid.species.all.cover.rain %>%
# arrange(desc(resid.cover.rain))#Of the species captured in both the seed rain and cover, make a plot showing which species were present in the seed rain in greater numbers than expected.
# resid.species.all.cover.rain %>%
# mutate(SPP6 = fct_reorder(SPP6, resid.cover.rain, .fun='median')) %>%
# ggplot(aes(x = SPP6, y = resid.cover.rain, color = Site)) +
# labs(x = "Species", y = "Residual", title = "Seed rain vs. Cover")+
# geom_point() +
# theme_classic() +
# theme(axis.text.x = element_text(angle = 70, hjust = 1))+
# geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1)+
# ylim(-5,5)+
# theme(axis.text = element_text( size = 12),
# axis.text.x = element_text( size = 5),
# axis.title = element_text( size = 16),
# strip.text = element_text(size = 16)) +
# annotate("text", x = 22, y = 3, label = "Overrepresented in seed rain")+annotate("text", x = 75, y = -4, label = "Underrepresented in seed rain")Here is a plot showing which species were present in the seed rain in greater numbers than expected based on local cover. Bars represent 95% confidence intervals. Absent bars for a species indicates it was only present in the cover and seed rain in one transect.
Remember that this figure only shows species that we were able to captured both in the cover and the seed rain at a transect. Does not show anything we caught only seeds or only cover at.
The figure does line up very well with the results of Maya’s seed predation study. Some of the species favored by rodents were Helianthus mollis, Sporobolus heterolepis, Desmodium canadense, Liatris aspera, and Baptisia bracteata. Based on the figure, these species or congeners also had far fewer seeds that you would expect based on local cover. Other species on this end of the spectrum are also well known to be favored by birds (e.g., American goldfinch and bobwhite). For example, Cirsium altissimum, Silphium integrifolium, Echinacea pallida, Strophostyles leiosperma, and Ambrosia artemisiifolia. Baptisia bracteata and Baptisia alba are also known to suffer heavy pre-dispersal predation from insects at Tucker Prairie (Haddock et al. 1982).
Lysimachia lanceolata is interesting since it has a tight relationship with a rare genera of bees, Macropis. It is unclear how dependent this species is on pollination by Macropis for reproduction.
Size, both plant and seed, are also likely important factors. Based on Moles et al 2004, small plants also tend to produce many small seeds (e.g., Juncus) while large plants (e.g., Baptisia) produce fewer, large seeds with tradeoffs existing between seed size and juvenile survival.
Primary dispersal mode is likely also important here - many species high in cover but low in captured seeds have fleshy fruits. For example, Potentilla simplex, Rhus glabra, Rhus copallinum, Rosa spp., and Fragraria virginica. I believe the only fleshy fruited species we captured was Rubus spp.
Life history strategy is also important, with annuals being grouped on the overrepresented end of the spectrum.
# resid.species.all.mean.cover.rain <- resid.species.all.cover.rain %>%
# group_by(SPP6) %>%
# summarize(mean.resid = mean(resid.cover.rain ),
# sd.resid = sd(resid.cover.rain),
# n.resid = n()) %>%
# mutate(se.resid = sd.resid / sqrt(n.resid),
# lower.ci.resid = mean.resid - 1.96*((sd.resid/sqrt(n.resid))),
# upper.ci.resid = mean.resid + 1.96*((sd.resid/sqrt(n.resid))))# resid.species.all.mean.cover.rain %>%
# mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
# ggplot(aes(x = SPP6, y = mean.resid)) +
# xlab("Species")+
# ylab("Residual") +
# geom_point() +
# geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
# theme_classic() +
# theme(axis.text.x = element_text(angle = 70, hjust = 1))+
# geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
# ylim(-2.5,2.5) +
# theme(axis.text = element_text( size = 10),
# axis.text.x = element_text( size = 5),
# axis.title = element_text( size = 16),
# strip.text = element_text(size = 16))+
# annotate("text", x = 30, y = 2.5, label = "Overrepresented in local seed rain")+
# annotate("text", x = 75, y = -2.5, label = "Underrepresented in local seed rain")
# ### Join the traits together
#
# resid.mean_cover_rain_traits <- left_join(resid.species.all.mean.cover.rain
# , traits)
#
# resid.mean_cover_rain_traits_reduced <- resid.mean_cover_rain_traits[, c("mean.resid", "SPP6", "Functional_Group", "Life_History_Lumped", "Mean_1_Seed_Mass_g", "Primary_Dispersal_Vector", "Provenance")]
#
#
# resid.mean_cover_rain_traits_reduced <- na.omit(resid.mean_cover_rain_traits_reduced)
#
# nrow(resid.mean_cover_rain_traits_reduced)# #### Life history vs. residuals
#
# ### Annual/biennials are more likely to be overrepresented
# ### Perennials are more likely to be underrepresented
#
# ggplot(data = resid.mean_cover_rain_traits_reduced, aes(x = Life_History_Lumped, y = mean.resid))+
# geom_boxplot()+
# geom_jitter()+
# geom_hline(yintercept = 0, color = "red")+
# labs(x = "Life History", y = "Pearson Residuals") +
# theme_classic() +
# theme(text=element_text(size=18), legend.key.size=unit(0.25, "cm"))
#
# #### Functional Group vs. residuals
#
# ### Forbs & Legumes don't look different from 0
# ### Graminoids are overrepresented
# ### Woody species are underrepresented
#
# ggplot(data = resid.mean_cover_rain_traits_reduced, aes(x = Functional_Group, y = mean.resid))+
# geom_boxplot()+
# geom_jitter()+
# geom_hline(yintercept = 0, color = "red")+
# labs(x = "Functional Group", y = "Pearson Residuals") +
# theme_classic() +
# theme(text=element_text(size=18), legend.key.size=unit(0.25, "cm"))+
# ylim(-2.1,2.1)
#
#
# #### Primary Dispersal Vector vs. residuals
#
# ## Animal and wind might be underrepresented
# ## Unassisted is neutral
#
# ggplot(data = resid.mean_cover_rain_traits_reduced, aes(x = Primary_Dispersal_Vector, y = mean.resid))+
# geom_boxplot()+
# geom_jitter()+
# geom_hline(yintercept = 0, color = "red")+
# labs(x = "Primary Dispersal Vector", y = "Pearson Residuals") +
# theme_classic() +
# theme(text=element_text(size=18), legend.key.size=unit(0.25, "cm"))+
# ylim(-2.1,2.1)
#
# #### Provenance vs. residuals
#
# ## Introduced has a slight trend toward over represented
#
# ggplot(data = resid.mean_cover_rain_traits_reduced, aes(x = Provenance, y = mean.resid))+
# geom_boxplot()+
# geom_jitter()+
# geom_hline(yintercept = 0, color = "red")+
# labs(x = "Provenance", y = "Pearson Residuals") +
# theme_classic() +
# theme(text=element_text(size=18), legend.key.size=unit(0.25, "cm"))+
# ylim(-2.1,2.1)
#
# #### Seed mass vs. residuals
#
# ## Introduced has a slight trend toward over represented
#
# ggplot(data = resid.mean_cover_rain_traits_reduced, aes(x = Mean_1_Seed_Mass_g, y = mean.resid))+
# scale_x_continuous(trans='log10') +
# geom_point()+
# geom_smooth(method = "lm") +
# theme_classic()+
# labs(x = "Seed Mass mg [log 10]", y = "Pearson Residuals")+
# theme_classic() +
# theme(text=element_text(size=18), legend.key.size=unit(0.25, "cm"))+
# geom_hline(yintercept = 0, color = "red")—-
Seed Bank and Vegetative Cover Exploratory Analysis
cover_bank_sum <- cover_bank %>%
group_by(SPP6) %>%
summarize(tot.cover = sum(Cover), tot.seeds = sum(Number_Seedlings))Germinated but didn’t observed in the cover
- Abutilon theophrasti
- Amaranthus tuberculatus
- Cardamine hirsuta
- Cardamine parviflora
- Digitaria sanguinalis
- Dipsacus laciniatus
- Euphorbia humistrata
- Hordeum pusillum
- Ipomoea hederacea
- Mollugo verticillata
- Panicum capillare
- Panicum dichotomiflorum
- Persicaria longiseta
- Platanus occidentalis (likely a GH contaminate)
- Plantago pusilla
- Poa chapmaniana
- Populus deltoides (likely a GH contaminate)
- Portulaca oleracea
- Setaria glauca
- Toxicodendron radicans
- Verbena hastata
- Veronica peregrina
- Verbascum thapsus
Almost everything we were able to germinate that was not in the local cover were weedy annual species. Some things like the Cardamine might have been in the cover but have such an early phenology that we missed them. The two tree species might be greenhouse contaminates since the roof was partially open and these trees were right next door.
Observed in the cover but didn’t germinate
Too many to list (112 taxa)
Clearly, only a small fraction of local species survive the transition from cover -> seed rain -> to germinable seed bank.
Few perennial species were present in the germinable seed bank.
Dissimilarity analysis
## # A tibble: 10,062 × 5
## # Groups: Site, Transect, SPP6 [7,644]
## Site Transect SPP6 Type Abundance
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 ABUTHE Seed_Bank 0
## 2 PFCA 1 1 ACAVIR Aboveground 11.5
## 3 PFCA 1 1 ACAVIR Seed_Bank 0
## 4 PFCA 1 1 ACERUB Aboveground 0
## 5 PFCA 1 1 ACHMIL Aboveground 0
## 6 PFCA 1 1 ACHMIL Seed_Bank 0
## 7 PFCA 1 1 AGAFAS Aboveground 0
## 8 PFCA 1 1 AGATEN Aboveground 0
## 9 PFCA 1 1 AGRGIG Aboveground 0
## 10 PFCA 1 1 AGRHYE Aboveground 0
## # ℹ 10,052 more rows# Code obtained and modified from Lauren and Anu's seed bank vs. vegetation paper
## Only included Bray-Curtis and Jaccard dissimilarity matrices since the others produced similar results to these two
##### run dissimilarity loops
cover_bank_dis_results <-matrix(nrow=0,ncol=4)
sites <-as.vector(unique(cover_bank_dis$Site))
for(s in 1:length(sites)){
temp <-subset(cover_bank_dis, Site==sites[s])
temp$SPP6 <- factor(as.character(temp$SPP6))
transects <-as.vector(unique(temp$Transect))
#for(b in 1:length(blocks)){
#temp_b<-subset(temp, block==blocks[b])
#plots<-as.vector(unique(temp_b$plot))
for(t in 1:length(transects)){
temp_t <- subset(temp, Transect==transects[t])
#wide form - with a row for above and belowground
plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
#relative abundances
trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"tot"))
names(trans_tot)[1]<-"trt"
trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
names(trans_pa)[1]<-"trt"
distance_bray<-vegdist(trans_tot[,-1], method = "bray")
distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")
new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1])
cover_bank_dis_results <-rbind( cover_bank_dis_results, new.row)
}
print(s/length(sites))
}## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1row.names(cover_bank_dis_results )<-NULL
cover_bank_dis_results <- data.frame(cover_bank_dis_results )
names(cover_bank_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard")##
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = cover_bank_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.214753 -0.058572 0.004864 0.067143 0.165558
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.826710 0.031169 26.524 < 2e-16 ***
## SitePFCA 1 -0.269145 0.042963 -6.265 3.47e-07 ***
## SitePFCA 2 -0.004166 0.042963 -0.097 0.923
## SitePFCA 3 0.045173 0.042963 1.051 0.300
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09351 on 35 degrees of freedom
## Multiple R-squared: 0.6665, Adjusted R-squared: 0.6379
## F-statistic: 23.31 on 3 and 35 DF, p-value: 1.801e-08## Anova Table (Type III tests)
##
## Response: dissimilarity_bray
## Sum Sq Df F value Pr(>F)
## (Intercept) 6.1511 1 703.508 < 2.2e-16 ***
## Site 0.6115 3 23.314 1.801e-08 ***
## Residuals 0.3060 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Cover_bank.mod.bray)
##
## $Site
## diff lwr upr p adj
## PFCA 1-TP -0.26914527 -0.38501257 -0.1532780 0.0000020
## PFCA 2-TP -0.00416572 -0.12003301 0.1117016 0.9996669
## PFCA 3-TP 0.04517336 -0.07069393 0.1610407 0.7206820
## PFCA 2-PFCA 1 0.26497955 0.15220261 0.3777565 0.0000016
## PFCA 3-PFCA 1 0.31431864 0.20154169 0.4270956 0.0000000
## PFCA 3-PFCA 2 0.04933908 -0.06343786 0.1621160 0.6433172##
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = cover_bank_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.108471 -0.051282 -0.007532 0.055739 0.133446
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.835304 0.022082 37.828 <2e-16 ***
## SitePFCA 1 -0.052772 0.030438 -1.734 0.0918 .
## SitePFCA 2 -0.006833 0.030438 -0.224 0.8237
## SitePFCA 3 0.016874 0.030438 0.554 0.5828
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06625 on 35 degrees of freedom
## Multiple R-squared: 0.1471, Adjusted R-squared: 0.07394
## F-statistic: 2.011 on 3 and 35 DF, p-value: 0.1302## Anova Table (Type III tests)
##
## Response: dissimilarity_jaccard
## Sum Sq Df F value Pr(>F)
## (Intercept) 6.2796 1 1430.9303 <2e-16 ***
## Site 0.0265 3 2.0114 0.1302
## Residuals 0.1536 35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = x)
##
## $Site
## diff lwr upr p adj
## PFCA 1-TP -0.052771996 -0.13485957 0.02931558 0.3220669
## PFCA 2-TP -0.006833257 -0.08892083 0.07525432 0.9959337
## PFCA 3-TP 0.016874135 -0.06521344 0.09896171 0.9447338
## PFCA 2-PFCA 1 0.045938739 -0.03395944 0.12583692 0.4192751
## PFCA 3-PFCA 1 0.069646131 -0.01025205 0.14954431 0.1059540
## PFCA 3-PFCA 2 0.023707392 -0.05619079 0.10360557 0.8538274Fig. - Bray-Curtis dissimilarity
## 95% confidence intervals
emmeans::emmeans(Cover_bank.mod.bray, "Site", type = "response")## Site emmean SE df lower.CL upper.CL
## TP 0.827 0.0312 35 0.763 0.890
## PFCA 1 0.558 0.0296 35 0.498 0.618
## PFCA 2 0.823 0.0296 35 0.763 0.883
## PFCA 3 0.872 0.0296 35 0.812 0.932
##
## Confidence level used: 0.95## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.558, 0.823, 0.872, 0.827)
lower <- c(0.498, 0.763, 0.812, 0.763)
upper <- c(0.618, 0.883, 0.932, 0.890)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Cover_bank_dissimilarity <- data.frame(site, response, lower, upper)
Fig. - Jaccard dissimilarity figure
## 95% confidence intervals
emmeans::emmeans(Cover_bank.mod.jac, "Site", type = "response")## Site emmean SE df lower.CL upper.CL
## TP 0.835 0.0221 35 0.790 0.880
## PFCA 1 0.783 0.0209 35 0.740 0.825
## PFCA 2 0.828 0.0209 35 0.786 0.871
## PFCA 3 0.852 0.0209 35 0.810 0.895
##
## Confidence level used: 0.95## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.783, 0.828, 0.852, 0.835)
lower <- c(0.740, .786, 0.810, 0.790)
upper <- c( 0.825, 0.871, 0.895, 0.880)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Cover_bank_jac_dissimilarity <- data.frame(site, response, lower, upper)
Immigrants
Proportion of new species in the seed bank not seen in local cover
cover_bank_dis_wide <- cover_bank_dis %>%
spread(SPP6, Abundance) %>%
mutate_all(~replace(., is.na(.), 0)) ## `mutate_all()` ignored the following grouping variables:
## • Columns `Site`, `Transect`
## ℹ Use `mutate_at(df, vars(-group_cols()), myoperation)` to silence the message.cover_bank_dis_long <- cover_bank_dis_wide %>%
gather(key = "SPP6", value = "Abundance", 4:199) %>%
spread(key = Type, Abundance)
# Get rid of columns that both have zeros
cover_bank_dis_long_unique <- cover_bank_dis_long[!(cover_bank_dis_long$Seed_Bank == 0 & cover_bank_dis_long$Aboveground == 0),]
# 0 = no seeds/seedlings and 1 = seeds/seedlings
cover_bank_dis_long_unique <- cover_bank_dis_long_unique %>%
mutate(Seed_Bank = if_else(Seed_Bank > 0, 1, 0)) %>%
mutate(Aboveground = if_else(Aboveground > 0, 1, 0))
# Make a new column that adds both the cover and seed bank presence/absence columns. Then remove rows that = 2 (shares both species)
cover_bank_dis_long_unique$shared <- cover_bank_dis_long_unique$Seed_Bank + cover_bank_dis_long_unique$Aboveground
cover_bank_dis_long_unique2 <- cover_bank_dis_long_unique[!(cover_bank_dis_long_unique$shared == 2),]
cover_bank_dis_long_unique2 <- cover_bank_dis_long_unique2[, -6]
# Count unique taxa per site and transect for the seed bank compared to the cover
Unique_species_bank_com_cover <- cover_bank_dis_long_unique2[, -4]
Unique_species_bank_com_cover <- Unique_species_bank_com_cover %>%
filter(Seed_Bank != 0 ) %>% group_by(Site, Transect) %>%
summarize(New_Bank = length(unique(SPP6)))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.# Count unique taxa per site and transect for the cover compared to the seed bank
Unique_species_cover_com_bank <- cover_bank_dis_long_unique2[, -5]
Unique_species_cover_com_bank <- Unique_species_cover_com_bank %>%
filter(Aboveground != 0 ) %>%
group_by(Site, Transect) %>%
summarize(New_Cover = length(unique(SPP6)))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.# Count the number of shared taxa between communities
cover_bank_dis_long_shared <- cover_bank_dis_long_unique %>%
filter(shared == 2) %>%
group_by(Site, Transect) %>%
summarize(Shared = length(unique(SPP6)))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.## Merge them together
Species_counts_cover_bank <- full_join(Unique_species_bank_com_cover, Unique_species_cover_com_bank) %>%
mutate_all(~replace(., is.na(.), 0)) ## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.Species_counts_cover_bank <- full_join(Species_counts_cover_bank, cover_bank_dis_long_shared) %>%
mutate_all(~replace(., is.na(.), 0)) ## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.Species_counts_cover_bank$Tot_Richness <- Species_counts_cover_bank$New_Bank + Species_counts_cover_bank$New_Cover + Species_counts_cover_bank$Shared
# Count the number of species in the seed rain for each transect
Species_counts_cover_bank$Tot_Cover <- Species_counts_cover_bank$New_Cover + Species_counts_cover_bank$Shared
# Count the number of species in the seed bank for each transect
Species_counts_cover_bank$Tot_Bank <- Species_counts_cover_bank$New_Bank + Species_counts_cover_bank$Shared
## Immigrants
Species_counts_cover_bank$Prop_Immigrant <- Species_counts_cover_bank$New_Bank / Species_counts_cover_bank$Tot_Cover### the oldest restoration has the greatest proportion of immigrants -> likely caused by a persistent seed bank formed from when this site was farmed prior to restoration. Weedy non-native species likely were displaced from the aboveground vegetation but remained present in the seed bank. This trend is not seen in Tucker Prairie which doesn't have the same past.
prop_new_species_cover_bank.plot <- ggplot() +
geom_jitter(data = Species_counts_cover_bank, aes(x = Site, y = Prop_Immigrant, color = Site), show.legend = FALSE, width = 0.2, size = 2.5) +
theme_classic()+
labs(x = "", y = "Prop. of new species in seed bank")+
scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"))+
ylim(-0.05,1)
prop_new_species_cover_bank.plot
Proportion of seeds in the seed bank not seen in local cover
names(cover_bank_dis_long) <- c("Site", "Transect", "SPP6", "Aboveground_Abundance", "Seed_Bank_Abundance")
cover_bank_dis_long_unique_abundance <- full_join(cover_bank_dis_long_unique, cover_bank_dis_long)## Joining with `by = join_by(Site, Transect, SPP6)`cover_bank_dis_long_shared_abundance <- cover_bank_dis_long_unique_abundance %>%
filter(shared == 2) %>%
group_by(Site, Transect) %>%
summarize(Shared_Seedlings = sum(Seed_Bank_Abundance))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.cover_bank_dis_long_immigrant_abundance <- cover_bank_dis_long_unique_abundance %>%
filter(shared == 1) %>%
group_by(Site, Transect) %>%
summarize(Immigrant_Seedlings = sum(Seed_Bank_Abundance))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.Prop_Immigrant_Seedlings <- full_join(cover_bank_dis_long_shared_abundance, cover_bank_dis_long_immigrant_abundance)## Joining with `by = join_by(Site, Transect)`Prop_Immigrant_Seedlings$Total_Seedlings <- Prop_Immigrant_Seedlings$Shared_Seedlings + Prop_Immigrant_Seedlings$Immigrant_Seedlings
Prop_Immigrant_Seedlings$Prop <- Prop_Immigrant_Seedlings$Immigrant_Seedlings / Prop_Immigrant_Seedlings$Total_Seedlingsprop_immigrant_seedlings_cover_bank.plot <- ggplot()+
geom_jitter(data=Prop_Immigrant_Seedlings, aes(x = Site, y = Prop, color = Site), show.legend = FALSE, width = .2, size = 2.5) +
theme_classic()+
labs(x = "", y = "Prop. of new seeds in seed bank", title = "Aboveground - Seed Bank")+
scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), plot.title = element_text(hjust = 0.5))+
ylim(-0.05,1)
prop_immigrant_seedlings_cover_bank.plot
Cover vs. Seed Bank Figures
By Site
# Drop Control for now
cover_bank <- cover_bank[!(cover_bank$Site%in%"CONTROL"),]ggplot(cover_bank, aes(x = Cover+1, y = Number_Seedlings+1))+
geom_point()+
scale_x_continuous(trans='log10') +
scale_y_continuous(trans='log10') +
facet_grid(col=vars(Site), labeller = labeller(Site = site.labs))+
geom_smooth(method = "lm") +
labs(x = "Veg. cover [Log 10 +1]", y = "No. Seedlings [Log 10 + 1]")+
theme(axis.text = element_text( size = 14),
axis.text.x = element_text( size = 16),
axis.title = element_text( size = 18),
legend.position = "none",
strip.text = element_text(size = 16))## `geom_smooth()` using formula = 'y ~ x'
Same Plot
ggplot(cover_bank, aes(x = Cover+1, y = Number_Seedlings+1, color = Site))+
geom_point()+
geom_smooth(method = "lm") +
scale_x_continuous(trans='log10') +
scale_y_continuous(trans='log10') +
labs(y = "Log (Total Germinated Seeds )", x = "Log (Total Cover )", title = "Seed bank vs cover")+
theme_classic() ## `geom_smooth()` using formula = 'y ~ x'
Residual plots
Fitting the Model
# # Remove everything that we didn't see both cover and seed bank for at a transect
#
# cover_bank_filtered <- cover_bank %>%
# filter(Cover > 0) %>%
# filter(Number_Seedlings > 0)
#
#
#
# #
# mod1 <- lmer(log(Number_Seedlings)~log(Cover)+(1|Site), data = cover_bank_filtered )
#
#
# plot(mod1)Again we get the streaks in the residual plots.
# p1<-ggplot(cover_bank_filtered , aes(x=log(Number_Seedlings))) +
# geom_histogram(color="black", fill="white")
# p1
#
# q1 <-ggplot(cover_bank_filtered , aes(x=log(Cover))) +
# geom_histogram(color="black", fill="white")
# q1 # resid.cover.bank <- residuals(mod1)
# #resid <- filtered.mod$residuals
#
#
# SPP6 <- cover_bank_filtered $SPP6
# Site <- cover_bank_filtered $Site
#
# resid.species.all.cover.bank <- data.frame(SPP6, resid.cover.bank, Site)
# resid.specieslall.cover.bank <- resid.species.all.cover.bank %>%
# arrange(desc(resid.cover.bank))#Of the species captured in both the seed bank and cover, make a plot showing which species were present in the seed bank in greater numbers than expected.
#
# resid.species.all.cover.bank %>%
# mutate(SPP6 = fct_reorder(SPP6, resid.cover.bank, .fun='median')) %>%
# ggplot(aes(x = SPP6, y = resid.cover.bank, color = Site)) +
# labs(x = "Species", y = "Residual", title = "Seed Bank vs. Cover")+
# geom_point() +
# theme_classic() +
# theme(axis.text.x = element_text(angle = 70, hjust = 1))+
# geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
# theme(axis.text = element_text( size = 10),
# axis.text.x = element_text( size = 8),
# axis.title = element_text( size = 16),
# strip.text = element_text(size = 16)) +
# ylim(-4,4)+
# annotate("text", x = 12, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 34, y = -4, label = "Underrepresented in seed bank")Looking at this residual plot, disturbance tolerant annual species tend to be over represented in the seed bank compared to the local cover. These species are most likely enriched in seed banks overtime and are persistent. For example, Digitaria ischaemum, Setaria faberi, Plantago virginica, and Juncus spp. Some over representation might be because of superb dispersal ability (Eupatorium serotinum and Agrostis hyemalis).
Life history is clearly important, with annuals being more present than perennials.
The seed bank and cover had the least amount of overlap with shared species.
# resid.species.all.mean.cover.bank <- resid.species.all.cover.bank %>%
# group_by(SPP6) %>%
# summarize(mean.resid = mean(resid.cover.bank),
# sd.resid = sd(resid.cover.bank),
# n.resid = n()) %>%
# mutate(se.resid = sd.resid / sqrt(n.resid),
# lower.ci.resid = mean.resid - 1.96*((sd.resid/sqrt(n.resid))),
# upper.ci.resid = mean.resid + 1.96*((sd.resid/sqrt(n.resid))))# resid.species.all.mean.cover.bank %>%
# mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
# ggplot(aes(x = SPP6, y = mean.resid)) +
# labs(x = "Species", y = "Residual", title = "Seed Bank vs. Cover")+
# geom_point() +
# geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
# theme_classic() +
# theme(axis.text.x = element_text(angle = 70, hjust = 1))+
# geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + theme(axis.text = element_text( size = 12),
# axis.text.x = element_text( size = 8),
# axis.title = element_text( size = 16),
# strip.text = element_text(size = 16))+
# ylim(-4,4)+
# annotate("text", x = 10, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 34, y = -4, label = "Underrepresented in seed bank")# Mass_coverage_cover.bank <- left_join(resid.species.all.mean.cover.bank, traits)
#
# Mass_coverage_cover.bank <- Mass_coverage_cover.bank %>%
# mutate(Seed_Mass = if_else(Mean_1_Seed_Mass_g > 0, 1, 0)) %>%
# mutate_all(~replace(., is.na(.), 0))# pane2 <- Mass_coverage_cover.bank %>%
# mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
# ggplot(aes(x = SPP6, y = mean.resid, color = as.factor(Seed_Mass))) +
# labs(x = "Species", y = "Residual", title = "Aboveground - Seed Bank")+
# ylab("Residual") +
# geom_point() +
# scale_color_manual(values = c("black", "blue"))+
# geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
# theme_classic() +
# theme(axis.text.x = element_text(angle = 70, hjust = 1))+
# geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
# ylim(-5,5) +
# theme(axis.text = element_text( size = 10),
# axis.text.x = element_text( size = 5),
# axis.title = element_text( size = 16),
# strip.text = element_text(size = 16))+
# annotate("text", x = 16, y = 3, label = "Overrepresented in local seed bank")+
# annotate("text", x = 32, y = -4, label = "Underrepresented in local seed bank")
#
# pane2# test2 <- lm(mean.resid~ log10(Mean_1_Seed_Mass_g*1000), Mass_coverage_cover.bank)
# summary(test2)
#
# ggplot(aes(x = Mean_1_Seed_Mass_g*1000, y = mean.resid), data = Mass_coverage_cover.bank)+
# geom_point()+
# geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1)+
# scale_x_continuous(trans='log10')+
# geom_smooth(method = "lm")+
# theme_classic()+
# labs(x = "Seed mass (mg) [log 10 scale]", y = "Residual")+
# theme(text=element_text(size=18), legend.key.size=unit(0.25, "cm"))—-
Seed Rain and Seed Bank Exploratory Analysis
rain_bank_sum <- rain_bank %>%
group_by(SPP6) %>%
summarize(tot.seeds = sum(Number_Seeds), tot.seedlings = sum(Number_Seedlings))Germinated but didn’t catch in the seed rain
- Abutilon theophrasti
- Digitalis sanguinalis
- Dipsacus laciniatus
- Euphorbia humistrata
- Lepidium densiflorum
- Panicum dichotomiflorum
- Plantago pusilla
- Poa chapmaniana
Mostly ruderal annual species.
Caught in the seed rain but didn’t germinate
Too many to list (61 taxa)
Dissimilarity analysis
## # A tibble: 10,062 × 5
## # Groups: Site, Transect, SPP6 [7,644]
## Site Transect SPP6 Type Abundance
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 ABUTHE Seed_Bank 0
## 2 PFCA 1 1 ACAVIR Aboveground 11.5
## 3 PFCA 1 1 ACAVIR Seed_Bank 0
## 4 PFCA 1 1 ACERUB Aboveground 0
## 5 PFCA 1 1 ACHMIL Aboveground 0
## 6 PFCA 1 1 ACHMIL Seed_Bank 0
## 7 PFCA 1 1 AGAFAS Aboveground 0
## 8 PFCA 1 1 AGATEN Aboveground 0
## 9 PFCA 1 1 AGRGIG Aboveground 0
## 10 PFCA 1 1 AGRHYE Aboveground 0
## # ℹ 10,052 more rows## # A tibble: 8,280 × 5
## # Groups: Site, Transect, SPP6 [5,640]
## Site Transect SPP6 Type Abundance
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 ABUTHE Seed_Bank 0
## 2 PFCA 1 1 ACAVIR Seed_Bank 0
## 3 PFCA 1 1 ACAVIR Seed_Rain 22
## 4 PFCA 1 1 ACHMIL Seed_Bank 0
## 5 PFCA 1 1 ACHMIL Seed_Rain 0
## 6 PFCA 1 1 AGASPP Seed_Rain 0
## 7 PFCA 1 1 AGRHYE Seed_Bank 0
## 8 PFCA 1 1 AGRHYE Seed_Rain 0
## 9 PFCA 1 1 ALOCAR Seed_Bank 0
## 10 PFCA 1 1 ALOCAR Seed_Rain 0
## # ℹ 8,270 more rows# Code obtained and modified from Lauren and Anu's seed bank vs. vegetation paper
## Only included Bray-Curtis and Jaccard dissimilarity matrices since the others produced similar results to these two
##### run dissimilarity loops
rain_bank_dis_results <-matrix(nrow=0,ncol=4)
sites <-as.vector(unique(rain_bank_dis$Site))
for(s in 1:length(sites)){
temp <-subset(rain_bank_dis, Site==sites[s])
temp$SPP6 <- factor(as.character(temp$SPP6))
transects <-as.vector(unique(temp$Transect))
for(t in 1:length(transects)){
temp_t <- subset(temp, Transect==transects[t])
#wide form - with a row for above and belowground
plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
#relative abundances
trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"tot"))
names(trans_tot)[1]<-"trt"
trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
names(trans_pa)[1]<-"trt"
distance_bray<-vegdist(trans_tot[,-1], method = "bray")
distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")
new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1])
rain_bank_dis_results <-rbind(rain_bank_dis_results, new.row)
}
print(s/length(sites))
}## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1row.names(rain_bank_dis_results)<-NULL
rain_bank_dis_results <- data.frame(rain_bank_dis_results)
names(rain_bank_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard")##
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = rain_bank_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30162 -0.06832 -0.00551 0.07355 0.26438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.68821 0.04061 16.948 <2e-16 ***
## SitePFCA 1 -0.13699 0.05743 -2.386 0.0224 *
## SitePFCA 2 0.10386 0.05743 1.809 0.0789 .
## SitePFCA 3 0.04302 0.05743 0.749 0.4586
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1284 on 36 degrees of freedom
## Multiple R-squared: 0.3458, Adjusted R-squared: 0.2913
## F-statistic: 6.344 on 3 and 36 DF, p-value: 0.001449## Anova Table (Type III tests)
##
## Response: dissimilarity_bray
## Sum Sq Df F value Pr(>F)
## (Intercept) 4.7363 1 287.2485 < 2.2e-16 ***
## Site 0.3138 3 6.3436 0.001449 **
## Residuals 0.5936 36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Rain_bank.mod.bray)
##
## $Site
## diff lwr upr p adj
## PFCA 1-TP -0.13699033 -0.29165040 0.01766974 0.0981225
## PFCA 2-TP 0.10385576 -0.05080431 0.25851583 0.2860455
## PFCA 3-TP 0.04302486 -0.11163521 0.19768493 0.8763881
## PFCA 2-PFCA 1 0.24084609 0.08618602 0.39550616 0.0009414
## PFCA 3-PFCA 1 0.18001519 0.02535512 0.33467526 0.0171829
## PFCA 3-PFCA 2 -0.06083090 -0.21549097 0.09382917 0.7160209##
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = rain_bank_dis_results)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.135047 -0.042571 0.002154 0.030196 0.128868
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.76505 0.02097 36.477 <2e-16 ***
## SitePFCA 1 -0.04026 0.02966 -1.357 0.183
## SitePFCA 2 -0.01470 0.02966 -0.496 0.623
## SitePFCA 3 0.04497 0.02966 1.516 0.138
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06632 on 36 degrees of freedom
## Multiple R-squared: 0.1949, Adjusted R-squared: 0.1278
## F-statistic: 2.905 on 3 and 36 DF, p-value: 0.04791## Anova Table (Type III tests)
##
## Response: dissimilarity_jaccard
## Sum Sq Df F value Pr(>F)
## (Intercept) 5.8529 1 1330.5892 < 2e-16 ***
## Site 0.0383 3 2.9051 0.04791 *
## Residuals 0.1584 36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Rain_bank.mod.jac)
##
## $Site
## diff lwr upr p adj
## PFCA 1-TP -0.04025516 -0.120137976 0.03962766 0.5336666
## PFCA 2-TP -0.01469815 -0.094580972 0.06518466 0.9595699
## PFCA 3-TP 0.04496722 -0.034915593 0.12485004 0.4387232
## PFCA 2-PFCA 1 0.02555700 -0.054325814 0.10543982 0.8244018
## PFCA 3-PFCA 1 0.08522238 0.005339566 0.16510520 0.0327730
## PFCA 3-PFCA 2 0.05966538 -0.020217438 0.13954820 0.2026732Fig. - Bray-Curtis dissimilarity
## 95% confidence intervals
emmeans::emmeans(Rain_bank.mod.bray, "Site", type = "response")## Site emmean SE df lower.CL upper.CL
## TP 0.688 0.0406 36 0.606 0.771
## PFCA 1 0.551 0.0406 36 0.469 0.634
## PFCA 2 0.792 0.0406 36 0.710 0.874
## PFCA 3 0.731 0.0406 36 0.649 0.814
##
## Confidence level used: 0.95## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.551, 0.792, 0.731, 0.688)
lower <- c(0.469, 0.710, 0.649, 0.606)
upper <- c(0.634, 0.874, 0.814, 0.771)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Rain_bank_dissimilarity <- data.frame(site, response, lower, upper)
Fig. - Jaccard dissimilarity figure
## 95% confidence intervals
emmeans::emmeans(Rain_bank.mod.jac, "Site", type = "response")## Site emmean SE df lower.CL upper.CL
## TP 0.765 0.021 36 0.723 0.808
## PFCA 1 0.725 0.021 36 0.682 0.767
## PFCA 2 0.750 0.021 36 0.708 0.793
## PFCA 3 0.810 0.021 36 0.767 0.853
##
## Confidence level used: 0.95## Make a new dataframe to hold the model estimates and 95% confidence interval bounds
response <- c(0.725, 0.750, 0.810, 0.765)
lower <- c(0.682, 0.708, 0.767, 0.723)
upper <- c(0.767, 0.793, 0.853, 0.808)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")
Rain_bank_jac_dissimilarity <- data.frame(site, response, lower, upper)
Immigrants
Proportion of new species in the seed bank not seen in the seed rain
rain_bank_dis_wide <- rain_bank_dis %>%
spread(SPP6, Abundance) %>%
mutate_all(~replace(., is.na(.), 0)) ## `mutate_all()` ignored the following grouping variables:
## • Columns `Site`, `Transect`
## ℹ Use `mutate_at(df, vars(-group_cols()), myoperation)` to silence the message.rain_bank_dis_long <- rain_bank_dis_wide %>%
gather(key = "SPP6", value = "Abundance", 4:144) %>%
spread(key = Type, Abundance)
# Get rid of columns that both have zeros
rain_bank_dis_long_unique <- rain_bank_dis_long[!(rain_bank_dis_long$Seed_Bank == 0 & rain_bank_dis_long$Seed_Rain == 0),]
# 0 = no seeds/seedlings and 1 = seeds/seedlings
rain_bank_dis_long_unique <- rain_bank_dis_long_unique %>%
mutate(Seed_Bank = if_else(Seed_Bank > 0, 1, 0)) %>%
mutate(Seed_Rain = if_else(Seed_Rain > 0, 1, 0))
# Make a new column that adds both the seed rain and seed bank presence/absence columns. Then remove rows that = 2 (shares both species)
rain_bank_dis_long_unique$shared <- rain_bank_dis_long_unique$Seed_Bank + rain_bank_dis_long_unique$Seed_Rain
rain_bank_dis_long_unique2 <- rain_bank_dis_long_unique[!(rain_bank_dis_long_unique$shared == 2),]
rain_bank_dis_long_unique2 <- rain_bank_dis_long_unique2[, -6]
# Count unique taxa per site and transect for the seed bank compared to the seed rain
Unique_species_bank_com_rain <- rain_bank_dis_long_unique2[, -5]
Unique_species_bank_com_rain <- Unique_species_bank_com_rain %>%
filter(Seed_Bank != 0 ) %>% group_by(Site, Transect) %>%
summarize(New_Bank = length(unique(SPP6)))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.# Count unique taxa per site and transect for the seed rain compared to the seed bank
Unique_species_rain_com_bank <- rain_bank_dis_long_unique2[, -4]
Unique_species_rain_com_bank <- Unique_species_rain_com_bank %>%
filter(Seed_Rain != 0 ) %>%
group_by(Site, Transect) %>%
summarize(New_Rain = length(unique(SPP6)))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.# Count the number of shared taxa between communities
rain_bank_dis_long_shared <- rain_bank_dis_long_unique %>%
filter(shared == 2) %>%
group_by(Site, Transect) %>%
summarize(Shared = length(unique(SPP6)))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.## Merge them together
Species_counts_rain_bank <- full_join(Unique_species_bank_com_rain, Unique_species_rain_com_bank) %>%
mutate_all(~replace(., is.na(.), 0)) ## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.Species_counts_rain_bank <- full_join(Species_counts_rain_bank, rain_bank_dis_long_shared) %>%
mutate_all(~replace(., is.na(.), 0)) ## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
## ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.Species_counts_rain_bank$Tot_Richness <- Species_counts_rain_bank$New_Bank + Species_counts_rain_bank$New_Rain + Species_counts_rain_bank$Shared
# Count the number of species in the seed rain for each transect
Species_counts_rain_bank$Tot_Rain <- Species_counts_rain_bank$New_Rain + Species_counts_rain_bank$Shared
# Count the number of species in the seed bank for each transect
Species_counts_rain_bank$Tot_Bank <- Species_counts_rain_bank$New_Bank + Species_counts_rain_bank$Shared
## Immigrants
Species_counts_rain_bank$Prop_Immigrant <- Species_counts_rain_bank$New_Bank / Species_counts_rain_bank$Tot_Rainprop_new_species_rain_bank.plot <- ggplot()+
geom_jitter(data=Species_counts_rain_bank, aes(x = Site, y = Prop_Immigrant, color = Site), show.legend = FALSE, width = .2, size = 2.5) +
theme_classic()+
labs(x = "", y = "Prop. of new species in seed bank")+
scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), plot.title = element_text(hjust = 0.5))+
ylim(-0.05,1)Proportion of seedlings in the seed bank not seen in seed rain
names(rain_bank_dis_long) <- c("Site", "Transect", "SPP6", "Seed_Bank_Abundance", "Seed_Rain_Abundance")
rain_bank_dis_long_unique_abundance <- full_join(rain_bank_dis_long_unique, rain_bank_dis_long)## Joining with `by = join_by(Site, Transect, SPP6)`rain_bank_dis_long_shared_abundance <- rain_bank_dis_long_unique_abundance %>%
filter(shared == 2) %>%
group_by(Site, Transect) %>%
summarize(Shared_Seedlings = sum(Seed_Bank_Abundance))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.rain_bank_dis_long_immigrant_abundance <- rain_bank_dis_long_unique_abundance %>%
filter(shared == 1) %>%
group_by(Site, Transect) %>%
summarize(Immigrant_Seedlings = sum(Seed_Bank_Abundance))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.Prop_Immigrant_Seedlings_Rain <- full_join(rain_bank_dis_long_shared_abundance, rain_bank_dis_long_immigrant_abundance)## Joining with `by = join_by(Site, Transect)`Prop_Immigrant_Seedlings_Rain$Total_Seedlings <- Prop_Immigrant_Seedlings_Rain$Shared_Seedlings + Prop_Immigrant_Seedlings_Rain$Immigrant_Seedlings
Prop_Immigrant_Seedlings_Rain$Prop <- Prop_Immigrant_Seedlings_Rain$Immigrant_Seedlings / Prop_Immigrant_Seedlings_Rain$Total_Seedlingsprop_immigrant_seedlings_rain_bank.plot <- ggplot()+
geom_jitter(data=Prop_Immigrant_Seedlings_Rain, aes(x = Site, y = Prop, color = Site), show.legend = FALSE, width = .2, size = 2.5) +
theme_classic()+
labs(x = "", y = "Prop. of new seeds in seed bank", title = "Seed Rain - Seed Bank")+
scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), plot.title = element_text(hjust = 0.5))+
ylim(-0.05,1)
prop_immigrant_seedlings_rain_bank.plot
Seed Rain vs. Seed Bank Figures
By Site
# Drop the control for now
rain_bank <- rain_bank[!(rain_bank$Site%in%"CONTROL"),]ggplot(rain_bank, aes(x = Number_Seeds+1, y = Number_Seedlings+1))+
geom_point() +
scale_x_continuous(trans='log10') +
scale_y_continuous(trans='log10') +
facet_grid(col=vars(Site), labeller = labeller(Site = site.labs))+
geom_smooth(method = "lm") +
labs(x = "No. Seeds [Log 10 + 1])", y = "No. Seedlings [Log 10 + 1]")+
theme(axis.text = element_text( size = 12),
axis.text.x = element_text( size = 12),
axis.title = element_text( size = 16),
legend.position = "none",
strip.text = element_text(size = 16))## `geom_smooth()` using formula = 'y ~ x'
Same Plot
ggplot(rain_bank, aes(x = Number_Seeds+1, y = Number_Seedlings+1, color = Site))+
geom_point()+
scale_x_continuous(trans='log10') +
scale_y_continuous(trans='log10') +
geom_smooth(method = "lm") +
labs(y = "Log (Total Captured Seed Rain )", x = "Log (Total Cover )", title = "Seed rain vs Seed Bank")+
theme_classic() ## `geom_smooth()` using formula = 'y ~ x'
Residual plots
#
# # Remove everything that we didn't see both the seed rain and bank at a transect
#
# rain_bank_filtered <- rain_bank %>%
# filter(Number_Seedlings> 0) %>%
# filter(Number_Seeds > 0)
#
# #
# mod3 <- lmer(log(Number_Seedlings)~log(Number_Seeds)+(1|Site), data = rain_bank_filtered)
#
#
# plot(mod3)Again we get the streaks in the residual vs. fitted values plot.
# p2<-ggplot(rain_bank_filtered , aes(x=log(Number_Seedlings))) +
# geom_histogram(color="black", fill="white")
# p2
#
# q2 <-ggplot(rain_bank_filtered , aes(x=log(Number_Seeds))) +
# geom_histogram(color="black", fill="white")
# q2
# #
# resid.rain.bank<- residuals(mod3)
#
# SPP6 <- rain_bank_filtered$SPP6
# Site <- rain_bank_filtered$Site
#
# resid.species.all.rain.bank <- data.frame(SPP6, resid.rain.bank, Site)
# resid.specieslall.rain.bank <- resid.species.all.rain.bank %>%
# arrange(desc(resid.rain.bank))#Of the species captured in both the seed rain and bank, make a plot showing which species were present in the seed bank in greater numbers than expected.
# resid.species.all.rain.bank %>%
# mutate(SPP6 = fct_reorder(SPP6, resid.rain.bank, .fun='median')) %>%
# ggplot(aes(x = SPP6, y = resid.rain.bank, color = Site)) +
# labs(x = "Species", y = "Residual", title = "Seed Rain vs. Seed Bank") +
# geom_point() +
# theme_classic() +
# theme(axis.text.x = element_text(angle = 70, hjust = 1))+
# geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1)+
# theme(axis.text = element_text( size = 12),
# axis.text.x = element_text( size = 8),
# axis.title = element_text( size = 16),
# strip.text = element_text(size = 16))+
# ylim(-4,4)+
# annotate("text", x = 12, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 40, y = -4, label = "Underrepresented in seed bank")# resid.species.all.mean.rain.bank <- resid.species.all.rain.bank %>%
# group_by(SPP6) %>%
# summarize(mean.resid = mean(resid.rain.bank),
# sd.resid = sd(resid.rain.bank),
# n.resid = n()) %>%
# mutate(se.resid = sd.resid / sqrt(n.resid),
# lower.ci.resid = mean.resid - 1.96*((sd.resid/sqrt(n.resid))),
# upper.ci.resid = mean.resid + 1.96*((sd.resid/sqrt(n.resid))))Interesting that some species underrepresented in the seed rain based on local cover are present in germinable seed bank (LYSLAN, CIRALT, STRLEI, EUTGYM, RATPIN). Does this suggest high mortality when transitionig into the seed bank and low mortality once in the seed bank?
Again, most species overrepresented in the seed bank based on local seed input are weedy and annual species. Most perennials fall on the underrepresented side of the spectrum.
Many of the native ones also seem to be species that do well in restorations (Lespedeza capitata, Ratibida pinnata, Andropogon gerardii, Sorghastrum nutans, Coreopsis lanceolata, Tradescantia ohiensis, Penstemon digitalis, Schizachyrium scoparium, Solidago altissima, Symphytrichum spp., etc.). It’s probably no coincidence that these species that recruit from seed rain also do well in restorations that are broadcast seeded. I wonder if restoration practices are accidentally selecting for species that recruit well from seed (limited dispersal opportunities and reinforcement from what initially recruits).
I think the most surprising species that germinated is Lysimachia lanceolata, which is not present in the PFCA restorations.
# resid.species.all.mean.rain.bank %>%
# mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
# ggplot(aes(x = SPP6, y = mean.resid)) +
# labs(x = "Species", y = "Residual", title = "Seed Rain vs. Seed Bank") +
# geom_point() +
# geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
# theme_classic() +
# theme(axis.text.x = element_text(angle = 70, hjust = 1))+
# geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + theme(axis.text = element_text( size = 12),
# axis.text.x = element_text( size = 8),
# axis.title = element_text( size = 16),
# strip.text = element_text(size = 16))+
# ylim(-4,4)+
# annotate("text", x = 12, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 40, y = -4, label = "Underrepresented in seed bank")# resid.species.all.mean.rain.bank %>%
# mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
# ggplot(aes(x = SPP6, y = mean.resid)) +
# labs(x = "Species", y = "Residual", title = "Seed Bank vs. Seed Rain")+
# geom_point() +
# geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
# theme_classic() +
# theme(axis.text.x = element_text(angle = 70, hjust = 1))+
# geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + theme(axis.text = element_text( size = 12),
# axis.text.x = element_text( size = 8),
# axis.title = element_text( size = 16),
# strip.text = element_text(size = 16))+
# ylim(-4,4)+
# annotate("text", x = 10, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 34, y = -4, label = "Underrepresented in seed bank")# Mass_coverage_rain.bank <- left_join(resid.species.all.mean.rain.bank, traits)
#
# Mass_coverage_rain.bank <- Mass_coverage_rain.bank %>%
# mutate(Seed_Mass = if_else(Mean_1_Seed_Mass_g > 0, 1, 0)) %>%
# mutate_all(~replace(., is.na(.), 0))# pane3 <- Mass_coverage_rain.bank %>%
# mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>%
# ggplot(aes(x = SPP6, y = mean.resid, color = as.factor(Seed_Mass))) +
# labs(x = "Species", y = "Residual", title = "Seed Rain - Seed Bank")+
# geom_point() +
# scale_color_manual(values = c("black", "blue"))+
# geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
# theme_classic() +
# theme(axis.text.x = element_text(angle = 70, hjust = 1))+
# geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
# ylim(-5,5) +
# theme(axis.text = element_text( size = 10),
# axis.text.x = element_text( size = 5),
# axis.title = element_text( size = 16),
# strip.text = element_text(size = 16),
# legend.position = "none")+
# annotate("text", x = 16, y = 3, label = "Overrepresented in local seed bank")+
# annotate("text", x = 32, y = -4, label = "Underrepresented in local seed bank")
#
# pane3# test3 <- lm(mean.resid~ log10(Mean_1_Seed_Mass_g*1000), Mass_coverage_rain.bank)
# summary(test3)
#
# ggplot(aes(x = Mean_1_Seed_Mass_g*1000, y = mean.resid), data = Mass_coverage_rain.bank)+
# geom_point()+
# geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1)+
# scale_x_continuous(trans='log10')+
# geom_smooth(method = "lm")+
# theme_classic()+
# labs(x = "Seed mass (mg) [log 10 scale]", y = "Residual")+
# theme(text=element_text(size=18), legend.key.size=unit(0.25, "cm"))—–
All Communities
### Merge everything
bank_to_merge_everything <- bank_to_merge_rain
for(i in 1:nrow(bank_to_merge_everything )){
# Group all Cerastium together
## CERGLO
if(bank_to_merge_everything [i,3] == "CERGLO"){bank_to_merge_everything [i,3] <- "CERSPP"}
# Group all Cyperus together
## DESSES
if(bank_to_merge_everything [i,3] == "CYPSTR"){bank_to_merge_everything [i,3] <- "CYPSPP"}
# Group all the Carex together
if(bank_to_merge_everything [i,3] == "CARFES"){bank_to_merge_everything [i,3] <- "CARSPP"}
# Group all Erigeron together
## ERIANN
## ERISTR
if(bank_to_merge_everything [i,3] == "ERIANN"){bank_to_merge_everything [i,3] <- "ERISPP"}
if(bank_to_merge_everything [i,3] == "ERISTR"){bank_to_merge_everything[i,3] <- "ERISPP"}
# Change Juncus interior to Juncus sp. (since we were able to ID JUNMAR in the cover but not JUNINT)
if(bank_to_merge_everything[i,3] == "JUNINT"){bank_to_merge_everything[i,3] <- "JUNSPP"}
}
bank_to_merge_everything <- bank_to_merge_everything[!(bank_to_merge_everything$Site%in%"CONTROL"),]
bank_to_merge_everything <- bank_to_merge_everything[!(bank_to_merge_everything$Site%in%"TP" & bank_to_merge_everything$Transect%in% "7"),]
bank_to_merge_everything$Type <- rep("Seed_Bank", nrow(bank_to_merge_everything))
names(bank_to_merge_everything) <- c("Site", "Transect", "SPP6", "Abundance", "Type")
cover_rain_everything <- cover_rain_dis
for(i in 1:nrow(cover_rain_everything )){
# Group all Cyperus together
## CYPSTR
## CYPECH
if(cover_rain_everything [i,3] == "CYPACU"){cover_rain_everything [i,3] <- "CYPSPP"}
if(cover_rain_everything [i,3] == "CYPSTR"){cover_rain_everything [i,3] <- "CYPSPP"}
if(cover_rain_everything [i,3] == "CYPECH"){cover_rain_everything [i,3] <- "CYPSPP"}
# Group all Oenothera together
## OENFIL
## OENBIE
if(cover_rain_everything [i,3] == "OENFIL"){cover_rain_everything [i,3] <- "OENSPP"}
if(cover_rain_everything [i,3] == "OENBIE"){cover_rain_everything [i,3] <- "OENSPP"}
}
everything_merged <- full_join(cover_rain_everything, bank_to_merge_everything)## Joining with `by = join_by(Site, Transect, SPP6, Type, Abundance)`unique(everything_merged$SPP6)## [1] "ACAVIR" "ACERUB" "ACHMIL" "AGASPP" "AGRGIG" "AGRHYE" "ALOCAR" "AMATUB"
## [9] "AMBART" "AMBPSI" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB"
## [17] "BAPAUS" "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD"
## [25] "CAPBUR" "CARDSP" "CARSPP" "CERSPP" "CHAFAS" "CHEALB" "CIRALT" "COMUMB"
## [33] "CONCAN" "CORLAN" "CORPAL" "CORSPP" "CORTRI" "CROSAG" "CYPSPP" "DALPUR"
## [41] "DAUCAR" "DESILL" "DESSPP" "DIAARM" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL"
## [49] "ELESPP" "ELYCAN" "ELYVIR" "ERASPE" "EREHIE" "ERISPP" "ERIVIL" "ERYYUC"
## [57] "EUPCOR" "EUPSPP" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA" "GALOBT"
## [65] "GENSPP" "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP" "HORPUS" "HYPHIR"
## [73] "HYPPUN" "IPOHED" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER" "LEPDEN"
## [81] "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR" "LEUVUL" "LIAPYC" "LIASPP"
## [89] "LIASQU" "LINSUL" "LOBSPI" "LYCAME" "LYSLAN" "LYTALA" "MEDLUP" "MELSPP"
## [97] "MOLVER" "MONFIS" "MUHGLA" "MUHSPP" "MYOVER" "OENSPP" "OXADIL" "PANCAP"
## [105] "PANVIR" "PARINT" "PENDIG" "PERLON" "PLAMAJ" "PLAOCC" "PLAVIR" "POAPRA"
## [113] "POLSPP" "POTSIM" "PYCPIL" "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA"
## [121] "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO"
## [129] "SCISPP" "SCLTRI" "SENMAR" "SETPAR" "SETSPP" "SILANT" "SILINT" "SILLAC"
## [137] "SISSPP" "SOLALT" "SOLCAR" "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT"
## [145] "SPHOBT" "SPILAC" "SPOCOM" "SPOHET" "STRLEI" "SYMNOV" "SYMORB" "SYMSPP"
## [153] "TAROFF" "THLARV" "TRAOHI" "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP"
## [161] "ULMPUM" "ULMSPP" "VERHAS" "VERHEL" "VEROSP" "VERSPP" "VIOSAG" "VIOSPP"
## [169] "VULOCT" "ZIZAUR" "CYPACU" "PLAPUS" "POACHA" "ABUTHE" "POPDEL" "EUPHUM"
## [177] "PANDIC" "POROLE" "VERTHA" "DIPLAC" "DIGSAN" "TOXRAD"### Sum everything by species by transect again
everything_merged <- everything_merged %>%
group_by(Site, Transect, SPP6, Type) %>%
summarize(Abundance = sum(Abundance))## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.everything_merged_wide <- everything_merged %>%
spread(key="SPP6", value="Abundance")
#Replace NA values with 0 for the cover
everything_merged_wide[is.na(everything_merged_wide)] <- 0
# Regroup for the cover
everything_merged_dis <- everything_merged_wide %>%
gather(key = "SPP6", value = "Abundance", 4:185)
## everything_merged_wide (community matrix, not transformed)
## everything_merged (long-form version of this community matrix)
everything_mat_labs <- everything_merged_wide[, c("Site", "Transect","Type")]
everything_mat <- everything_merged_wide[, 4:185]
everything_mat_hell <- decostand(everything_mat, "hellinger")
everything_mat_hell_bray <- vegdist(everything_mat_hell, method = "bray")PCoA
library(ecodist)## Registered S3 method overwritten by 'ecodist':
## method from
## dim.dist proxy##
## Attaching package: 'ecodist'## The following object is masked from 'package:vegan':
##
## mantelbray_curtis_pcoa <- ecodist::pco(everything_mat_hell_bray )
# Creates a data frame from principal coordinates
bray_pcoa_df <- data.frame(pcoa1 = bray_curtis_pcoa$vectors[,1],
pcoa2 = bray_curtis_pcoa$vectors[,2])
labs_bray_pcoa_df <- data.frame(everything_mat_labs, bray_pcoa_df)
labs_bray_pcoa_df_centroid <- labs_bray_pcoa_df %>%
group_by(Site, Type) %>%
summarize(centroid.1 = mean(pcoa1), centroid.2 = mean(pcoa2))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.# Creates the plot
bray_plot <- ggplot(data = labs_bray_pcoa_df, aes(x=pcoa1, y=pcoa2, color = Site, shape = Type)) +
geom_point(alpha = .7, cex = 2) +
geom_point(data = labs_bray_pcoa_df_centroid , aes(x = centroid.1, y = centroid.2), cex = 4)+
labs(x = "PC1",
y = "PC2",
title = "Bray-curtis PCoA with Hellinger transformation") +
scale_shape_manual(name = "Type", breaks = c("Aboveground", "Seed_Bank", "Seed_Rain"),
labels = c("Aboveground", "Seed Bank", "Seed Rain"),
values = c(16, 17, 15))+
scale_color_manual(name = "Site", breaks = c("PFCA 1", "PFCA 2", "PFCA 3", "TP"),
labels = c("Young", "Middle", "Old", "Remnant"),
values = c("#E69F00", "#009E73", "#56B4E9", "#0072B2"))+
theme(title = element_text(size = 12)) +
theme_classic()
bray_plot
test <- adonis2(everything_mat_hell_bray ~ Site*Type, method= "bray", data = everything_mat_labs)
summary(test)## Df SumOfSqs R2 F
## Min. : 2.0 Min. : 3.778 Min. :0.1093 Min. : 4.62
## 1st Qu.: 3.0 1st Qu.: 5.656 1st Qu.:0.1636 1st Qu.:12.69
## Median : 6.0 Median :10.837 Median :0.3134 Median :20.75
## Mean : 46.4 Mean :13.833 Mean :0.4000 Mean :17.29
## 3rd Qu.:105.0 3rd Qu.:14.310 3rd Qu.:0.4138 3rd Qu.:23.63
## Max. :116.0 Max. :34.581 Max. :1.0000 Max. :26.51
## NA's :2
## Pr(>F)
## Min. :0.001
## 1st Qu.:0.001
## Median :0.001
## Mean :0.001
## 3rd Qu.:0.001
## Max. :0.001
## NA's :2library(pairwiseAdonis)## Loading required package: cluster### Ok good everything is different from everything
### there are a lot more similarities between TP and PFCA 3 because we had to lump a lot of things to genus level due to discrepencies in being able to ID everything
pairwise.adonis2(everything_mat_hell_bray ~ Site*Type, method= "bray", data = everything_mat_labs)## $parent_call
## [1] "everything_mat_hell_bray ~ Site * Type , strata = Null , permutations 999"
##
## $`PFCA 1_vs_PFCA 2`
## Df SumOfSqs R2 F Pr(>F)
## Site 1 1.9877 0.15853 16.8705 0.001 ***
## Type 2 3.4434 0.27462 14.6126 0.001 ***
## Site:Type 2 0.7451 0.05943 3.1621 0.001 ***
## Residual 54 6.3625 0.50742
## Total 59 12.5388 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`PFCA 1_vs_PFCA 3`
## Df SumOfSqs R2 F Pr(>F)
## Site 1 5.1443 0.30141 36.2324 0.001 ***
## Type 2 2.8476 0.16685 10.0281 0.001 ***
## Site:Type 2 1.4084 0.08252 4.9597 0.001 ***
## Residual 54 7.6670 0.44922
## Total 59 17.0673 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`PFCA 1_vs_TP`
## Df SumOfSqs R2 F Pr(>F)
## Site 1 6.0660 0.37271 48.6763 0.001 ***
## Type 2 2.3224 0.14269 9.3179 0.001 ***
## Site:Type 2 1.5316 0.09411 6.1452 0.001 ***
## Residual 51 6.3556 0.39050
## Total 56 16.2757 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`PFCA 2_vs_PFCA 3`
## Df SumOfSqs R2 F Pr(>F)
## Site 1 2.6222 0.16230 17.8013 0.001 ***
## Type 2 4.4245 0.27385 15.0182 0.001 ***
## Site:Type 2 1.1556 0.07152 3.9225 0.001 ***
## Residual 54 7.9544 0.49233
## Total 59 16.1567 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`PFCA 2_vs_TP`
## Df SumOfSqs R2 F Pr(>F)
## Site 1 3.8134 0.24391 29.2764 0.001 ***
## Type 2 3.6334 0.23240 13.9473 0.001 ***
## Site:Type 2 1.5447 0.09880 5.9296 0.001 ***
## Residual 51 6.6430 0.42489
## Total 56 15.6345 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`PFCA 3_vs_TP`
## Df SumOfSqs R2 F Pr(>F)
## Site 1 2.0978 0.13719 13.4620 0.001 ***
## Type 2 4.0498 0.26485 12.9940 0.001 ***
## Site:Type 2 1.1957 0.07820 3.8365 0.001 ***
## Residual 51 7.9476 0.51976
## Total 56 15.2909 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## attr(,"class")
## [1] "pwadstrata" "list"Trajectory Analysis
library(ecotraj)## Loading required package: RcppD = everything_mat_hell_bray
labs_bray_pcoa_df$Transect = factor(labs_bray_pcoa_df$Transect, levels=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
labs_bray_pcoa_df <- labs_bray_pcoa_df %>%
unite("Site_Transect", c("Site", "Transect"), na.rm = TRUE, remove = FALSE)
sites <- factor(labs_bray_pcoa_df$Site_Transect)
sites <- as.numeric(sites)
labs_bray_pcoa_df$Type <- factor(labs_bray_pcoa_df$Type, levels = c("Aboveground", "Seed_Rain", "Seed_Bank"))
surveys <- as.numeric(labs_bray_pcoa_df$Type)
### you are going to have to add site to the community level label PFCA1_Aboveground, PFCA1_Seed_Rain, PFCA1_Seed_Bank, ... to get 12 categories.
### Also I believe transect should be site and survey should be site/community type
# 1 = Aboveground, 2 = Seed Rain, 3 = Seed Bank
Young <- rep("#E69F00", 10)
Middle <- rep("#009E73", 10)
Old <- rep( "#56B4E9", 10)
Remnant <- rep("#0072B2", 10)
trajectoryPCoA(D, sites, traj.colors=c(Young, Middle, Old, Remnant), surveys, lwd = 2,
survey.labels = T)
legend(-1.1, .5, col=c("#E69F00", "#009E73", "#56B4E9","#0072B2"),
legend=c("Young", "Middle", "Old", "Remnant"), bty="n", lty=1, lwd = 2)
Site_labels <- c(rep("Young", 10), rep("Middle", 10), rep("Old", 10), rep("Remnant", 9))
seg.dist <- trajectoryLengths(D, sites, surveys)
Directionality <- trajectoryDirectionality(D, sites, surveys)
trajectory.stats <- cbind(Site_labels, seg.dist, Directionality)
trajectory.stats2 <- trajectory.stats %>%
gather(key = "Segment", value = "Length", c("S1", "S2", "Trajectory"))trajectory.stats2$Site_labels = factor(trajectory.stats2$Site_labels, levels=c("Young", "Middle", "Old", "Remnant"))
mean.distance <- trajectory.stats2 %>%
group_by(Site_labels, Segment) %>%
summarize(mean.length = mean(Length), sd.length = sd(Length))## `summarise()` has grouped output by 'Site_labels'. You can override using the
## `.groups` argument.ggplot(trajectory.stats2 , aes(x = Segment, y = Length, fill = Site_labels)) +
geom_boxplot(show.legend = FALSE) +
theme_classic() +
labs(x = "Transition",
y = "Distance") +
scale_x_discrete(name = "Type", breaks = c("S1", "S2", "Trajectory"),
labels = c("Above to SR", "SR to SB", "Above to SB"))+
scale_fill_manual(name = "Site", breaks = c("Young", "Middle", "Old", "Remnant"),
labels = c("Young", "Middle", "Old", "Remnant"),
values = c("#E69F00", "#009E73", "#56B4E9", "#0072B2"))+
theme(text=element_text(size=18), legend.key.size=unit(1, "cm"))+
theme_classic()
ggplot(trajectory.stats2 , aes(x = Site_labels, y = Directionality, fill = Site_labels)) +
labs(x = "Site", y = "Directionality")+
geom_boxplot(show.legend = FALSE) +
theme_classic() +
scale_fill_manual(name = "Site", breaks = c("Young", "Middle", "Old", "Remnant"),
labels = c("Young", "Middle", "Old", "Remnant"),
values = c("#E69F00", "#009E73", "#56B4E9", "#0072B2"))+
theme(text=element_text(size=18), legend.key.size=unit(1, "cm"))+
theme_classic()
Ternary Plot
## Step 1: Convert everything into proportions so everything is on the same scale
everything_merged_dis # <- dataset that contains every community combined AND with zeros for species that were not observered at a transect for each community## # A tibble: 21,294 × 5
## # Groups: Site, Transect [39]
## Site Transect Type SPP6 Abundance
## <chr> <fct> <chr> <chr> <dbl>
## 1 PFCA 1 1 Aboveground ABUTHE 0
## 2 PFCA 1 1 Seed_Bank ABUTHE 0
## 3 PFCA 1 1 Seed_Rain ABUTHE 0
## 4 PFCA 1 10 Aboveground ABUTHE 0
## 5 PFCA 1 10 Seed_Bank ABUTHE 0
## 6 PFCA 1 10 Seed_Rain ABUTHE 0
## 7 PFCA 1 2 Aboveground ABUTHE 0
## 8 PFCA 1 2 Seed_Bank ABUTHE 0
## 9 PFCA 1 2 Seed_Rain ABUTHE 0
## 10 PFCA 1 3 Aboveground ABUTHE 0
## # ℹ 21,284 more rows### Step 1.1 - make a dataset containing the totals for each community type and species
community_type_totals <- everything_merged_dis %>%
group_by(Site, Type, SPP6) %>%
summarize(tot.abundance = sum(Abundance))## `summarise()` has grouped output by 'Site', 'Type'. You can override using the
## `.groups` argument.community_type_lumped <- everything_merged_dis %>%
group_by(Site, Type) %>%
summarize(everything = sum(Abundance))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.everything_merged_dis_totals <- left_join(community_type_totals, community_type_lumped)## Joining with `by = join_by(Site, Type)`everything_merged_dis_totals$percent <- (everything_merged_dis_totals$tot.abundance / everything_merged_dis_totals$everything)*100### Note that ggtern does the normalizing for you!
library(ggtern)## Registered S3 methods overwritten by 'ggtern':
## method from
## grid.draw.ggplot ggplot2
## plot.ggplot ggplot2
## print.ggplot ggplot2## --
## Remember to cite, run citation(package = 'ggtern') for further info.
## --##
## Attaching package: 'ggtern'## The following objects are masked from 'package:ggplot2':
##
## aes, annotate, ggplot, ggplot_build, ggplot_gtable, ggplotGrob,
## ggsave, layer_data, theme_bw, theme_classic, theme_dark,
## theme_gray, theme_light, theme_linedraw, theme_minimal, theme_voidlibrary(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:mosaic':
##
## do## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(readr)
everything_merged_dis_totals_test <- everything_merged_dis_totals[, c("Site", "Type", "SPP6", "tot.abundance")]
everything_merged_dis_totals_wide <- everything_merged_dis_totals_test %>%
spread(key="Type", value="tot.abundance")
#Replace NA values with 0 for the cover
tern.plot.data <- everything_merged_dis_totals_wide %>%
filter(!(Seed_Rain == 0 & Aboveground == 0 & Seed_Bank == 0))
### Tucker Prairie
tern.plot.data_TP <- tern.plot.data %>%
filter(Site == "TP")
DATA <- data.frame(x = c(1,0,0),
y = c(0,1,0),
z = c(0,0,1),
xend = c(0,.5,.5),
yend = c(.5,0,.5),
zend = c(.5,.5,0),
Series = c("yz","xz","xy"))
ggtern(data=tern.plot.data_TP,aes(x=Aboveground,y=Seed_Bank, z=Seed_Rain))+
geom_mask() +
geom_segment(data=DATA,aes(x,y,z,xend=xend,yend=yend,zend=zend, color = Series), size = .75, linetype = "dotted", show.legend = FALSE)+
scale_color_manual(values=c("darkgreen","darkblue","darkred"))+
# geom_point(position= position_jitter_tern(x=0.02, y=0.02, z=0.02))+
geom_text(data = tern.plot.data_TP, aes(label = SPP6), size = 2)+
theme_rgbw()+
labs(title = "Tucker Prairie")## Warning in geom_segment(data = DATA, aes(x, y, z, xend = xend, yend = yend, :
## Ignoring unknown aesthetics: z and zend
#### Prairie Fork 1
tern.plot.data_PFCA_1 <- tern.plot.data %>%
filter(Site == "PFCA 1")
ggtern(data=tern.plot.data_PFCA_1,aes(x=Aboveground,y=Seed_Bank, z=Seed_Rain))+
geom_mask() +
geom_segment(data=DATA,aes(x,y,z,xend=xend,yend=yend,zend=zend, color = Series), size = .75, linetype = "dotted", show.legend = FALSE)+
scale_color_manual(values=c("darkgreen","darkblue","darkred"))+
#geom_point(position= position_jitter_tern(x=0.02, y=0.02, z=0.02))+
geom_text(data = tern.plot.data_PFCA_1, aes(label = SPP6), size = 2)+
theme_rgbw()+
labs(title = "Young")## Warning in geom_segment(data = DATA, aes(x, y, z, xend = xend, yend = yend, :
## Ignoring unknown aesthetics: z and zend
#### Prairie Fork 2
tern.plot.data_PFCA_2 <- tern.plot.data %>%
filter(Site == "PFCA 2")
ggtern(data=tern.plot.data_PFCA_2,aes(x=Aboveground,y=Seed_Bank, z=Seed_Rain))+
geom_mask() +
geom_segment(data=DATA,aes(x,y,z,xend=xend,yend=yend,zend=zend, color = Series), size = .75, linetype = "dotted", show.legend = FALSE)+
scale_color_manual(values=c("darkgreen","darkblue","darkred"))+
# geom_point(position= position_jitter_tern(x=0.02, y=0.02, z=0.02))+
geom_text(data = tern.plot.data_PFCA_2, aes(label = SPP6), size = 2)+
theme_rgbw()+
labs(title = "Middle")## Warning in geom_segment(data = DATA, aes(x, y, z, xend = xend, yend = yend, :
## Ignoring unknown aesthetics: z and zend
#### Prairie Fork 3
tern.plot.data_PFCA_3 <- tern.plot.data %>%
filter(Site == "PFCA 3")
ggtern(data=tern.plot.data_PFCA_3,aes(x=Aboveground,y=Seed_Bank, z=Seed_Rain))+
geom_mask() +
geom_segment(data=DATA,aes(x,y,z,xend=xend,yend=yend,zend=zend, color = Series), size = .75, linetype = "dotted", show.legend = FALSE)+
scale_color_manual(values=c("darkgreen","darkblue","darkred"))+
#geom_point(position= position_jitter_tern(x=0.02, y=0.02, z=0.02))+
geom_text(data = tern.plot.data_PFCA_3, aes(label = SPP6), size = 2)+
theme_rgbw()+
labs(title = "Old")## Warning in geom_segment(data = DATA, aes(x, y, z, xend = xend, yend = yend, :
## Ignoring unknown aesthetics: z and zend
—- Multipaned figures
# Dissimiliarity
plot_grid(bray_curtis_cover_rain.plot,bray_curtis_rain_bank.plot, bray_curtis_cover_bank.plot, ncol = 3)
# Sig differences for all except for SB-Veg Jaccard
## Still need to do pairwise differences (contrasts? Tukey HSD?)plot_grid(prop_immigrant_seeds_cover_rain.plot, prop_immigrant_seedlings_rain_bank.plot,prop_immigrant_seedlings_cover_bank.plot, ncol = 3)
Just_seedlings <- everything_merged %>%
filter(Type == "Seed_Bank") %>%
group_by(Site, Transect) %>%
summarize(tot.seedlings = sum(Abundance))## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.Just_seedlings_mean <- Just_seedlings %>%
group_by(Site) %>%
summarize(mean.seedlings = mean(tot.seedlings), sd.seedlings = sd(tot.seedlings))
ggplot(Just_seedlings, aes(x= Site, y = tot.seedlings, color = Site))+
geom_jitter(width =0.2, size = 2, show.legend = FALSE)+
geom_point(data = Just_seedlings_mean, aes(x = Site, y = mean.seedlings),size = 3.5, color = "black")+
geom_errorbar(data = Just_seedlings_mean, aes(x= Site, y = mean.seedlings, ymin = mean.seedlings - sd.seedlings, ymax = mean.seedlings + sd.seedlings), width = 0.15, cex = .75, color = "black")+
theme_classic()+
labs(x= "Site", y = "No. Seedlings")+
theme(text=element_text(size=18), legend.key.size=unit(1, "cm"), plot.title = element_text(hjust = 0.5))+
scale_color_manual(name = "Site", breaks = c("PFCA 1", "PFCA 2", "PFCA 3", "TP"),
labels = c("Young", "Middle", "Old", "Remnant"),
values = c("#1f78b4", "#a6cee3", "#b2df8a","#33a02c"))+
scale_x_discrete(name = "Site", breaks = c("PFCA 1", "PFCA 2", "PFCA 3", "TP"),
labels = c("Young", "Middle", "Old", "Remnant"))